source: readTracFiles_mod.f90 @ 22

Last change on this file since 22 was 22, checked in by dcugnet, 2 years ago

Fix needed for gfortran 4.8.5 compiler (JYG)

File size: 100.4 KB
Line 
1MODULE readTracFiles_mod
2
3  USE strings_mod,    ONLY: msg, testFile,  strFind, strStack, strReduce,  strHead, strCount, find, fmsg, reduceExpr, &
4             removeComment, cat, checkList, str2int, strParse, strReplace, strTail, strIdx, maxlen, test, dispTable, get_in
5  USE trac_types_mod, ONLY: trac_type, isot_type, keys_type
6
7  IMPLICIT NONE
8
9  PRIVATE
10
11  PUBLIC :: initIsotopes, maxlen, trac_type, isot_type, keys_type
12  PUBLIC :: readTracersFiles, indexUpdate, setGeneration             !--- TOOLS ASSOCIATED TO TRACERS  DESCRIPTORS
13  PUBLIC :: readIsotopesFile                                         !--- TOOLS ASSOCIATED TO ISOTOPES DESCRIPTORS
14  PUBLIC :: getKey_init, getKey, fGetKey, setDirectKeys              !--- GET/SET KEYS FROM/TO tracers & isotopes
15
16  PUBLIC :: addPhase, new2oldName,  getPhase, &                      !--- FUNCTIONS RELATED TO THE PHASES
17            delPhase, old2newName, getiPhase, &                      !--- + ASSOCIATED VARIABLES
18            known_phases, old_phases, phases_sep, phases_names, nphases
19
20  PUBLIC :: oldH2OIso, newH2OIso                                     !--- NEEDED FOR BACKWARD COMPATIBILITY (OLD traceur.def)
21
22  PUBLIC :: tran0, idxAncestor, ancestor                             !--- GENERATION 0 TRACER + TOOLS FOR GENERATIONS
23  PUBLIC :: maxTableWidth
24!------------------------------------------------------------------------------------------------------------------------------
25  TYPE :: dataBase_type                                              !=== TYPE FOR TRACERS SECTION
26    CHARACTER(LEN=maxlen)  :: name                                   !--- Section name
27    TYPE(trac_type), ALLOCATABLE :: trac(:)                          !--- Tracers descriptors
28  END TYPE dataBase_type
29!------------------------------------------------------------------------------------------------------------------------------
30  INTERFACE getKey
31    MODULE PROCEDURE getKeyByName_s1, getKeyByName_i1, getKeyByName_r1, getKeyByName_sm, getKeyByName_im, getKeyByName_rm
32  END INTERFACE getKey
33!------------------------------------------------------------------------------------------------------------------------------
34  INTERFACE fGetKey;       MODULE PROCEDURE fgetKeyByIndex_s1, fgetKeyByName_s1; END INTERFACE fGetKey
35  INTERFACE tracersSubset; MODULE PROCEDURE trSubset_Indx, trSubset_Name, trSubset_gen0Name; END INTERFACE tracersSubset
36  INTERFACE idxAncestor;   MODULE PROCEDURE idxAncestor_1, idxAncestor_m; END INTERFACE idxAncestor
37  INTERFACE    ancestor;   MODULE PROCEDURE    ancestor_1,    ancestor_m; END INTERFACE    ancestor
38  INTERFACE    addPhase;   MODULE PROCEDURE   addPhase_s1,   addPhase_sm,   addPhase_i1,   addPhase_im; END INTERFACE addPhase
39  INTERFACE old2newName;   MODULE PROCEDURE old2newName_1, old2newName_m; END INTERFACE old2newName
40  INTERFACE new2oldName;   MODULE PROCEDURE new2oldName_1, new2oldName_m; END INTERFACE new2oldName
41!------------------------------------------------------------------------------------------------------------------------------
42
43  !=== MAIN DATABASE: files sections descriptors
44  TYPE(dataBase_type), SAVE, ALLOCATABLE, TARGET :: dBase(:)
45
46  !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
47  CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'           !--- Default transporting fluid
48  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlir'          !--- Old phases for water (no separator)
49  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsr'          !--- Known phases initials
50  INTEGER,               PARAMETER :: nphases=LEN_TRIM(known_phases) !--- Number of phases
51  CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &        !--- Known phases names
52                                = ['gaseous', 'liquid ', 'solid  ', 'cloud  ']
53  CHARACTER(LEN=1), SAVE :: phases_sep  =  '_'                       !--- Phase separator
54  LOGICAL,          SAVE :: tracs_merge = .TRUE.                     !--- Merge/stack tracers lists
55  LOGICAL,          SAVE :: lSortByGen  = .TRUE.                     !--- Sort by growing generation
56
57  !--- KEPT JUST TO MANAGE OLD WATER ISOTOPES NAMES
58  !--- Apart from that context, on limitaion on isotopes names (as long as they have a corresponding line in isotopes_params.def)
59  CHARACTER(LEN=maxlen), SAVE :: oldH2OIso(5) = ['eau',     'HDO',     'O18',     'O17',     'HTO'    ]
60  CHARACTER(LEN=maxlen), SAVE :: newH2OIso(5) = ['H2[16]O', 'H[2]HO ', 'H2[18]O', 'H2[17]O', 'H[3]HO ']
61
62  !=== LOCAL COPIES OF THE TRACERS AND ISOTOPES DESCRIPTORS, USED BY getKey (INITIALIZED IN getKey_init)
63  TYPE(trac_type), ALLOCATABLE, TARGET, SAVE ::  tracers(:)
64  TYPE(isot_type), ALLOCATABLE, TARGET, SAVE :: isotopes(:)
65
66  INTEGER,    PARAMETER :: maxTableWidth = 192                       !--- Maximum width of a table displayed with "dispTable"
67  CHARACTER(LEN=maxlen) :: modname
68
69CONTAINS
70
71!==============================================================================================================================
72!==============================================================================================================================
73!=== READ ONE OR SEVERAL TRACER FILES AND FILL A "tr" TRACERS DESCRIPTOR DERIVED TYPE.
74!=== THE RETURN VALUE fType DEPENDS ON WHAT IS FOUND:
75!===  0: NO ADEQUATE FILE FOUND ; DEFAULT VALUES MUST BE USED
76!===  1: AN "OLD STYLE" TRACERS FILE "traceur.def":
77!===    First line: <nb tracers>     Other lines: <hadv> <vadv> <tracer name> [<parent name>]
78!===  2: A  "NEW STYLE" TRACERS FILE  "tracer.def" WITH SEVERAL SECTIONS.
79!===  3: SEVERAL  "  "  TRACERS FILES "tracer_<component>.def" WITH A SINGLE SECTION IN EACH.
80!=== REMARKS:
81!===  * EACH SECTION BEGINS WITH A "&<section name> LINE
82!===  * DEFAULT VALUES FOR ALL THE SECTIONS OF THE FILE ARE DEFINED IN THE SPECIAL SECTION "&default"
83!===  * EACH SECTION LINE HAS THE STRUCTURE:  <name(s)>  <key1>=<value1> <key2>=<value2> ...
84!===  * SO FAR, THE DEFINED KEYS ARE: parent, phases, hadv, vadv, type
85!===  * <name> AND <parent> CAN BE LISTS OF COMA-SEPARATED TRACERS ; THE ROUTINE EXPAND THESE FACTORIZATIONS.
86!=== FUNCTION RETURN VALUE "lerr" IS FALSE IN CASE SOMETHING WENT WRONG.
87!=== ABOUT THE KEYS:
88!     * The "keys" component (of type keys_type) is in principle enough to store everything we could need.
89!     But some variables are stored as direct-access keys to make the code more readable and because they are used often.
90!     * Most of the direct-access keys are set in this module, but some are not (longName, iadv, isAdvected for now).
91!     * Some of the direct-access keys must be updated (using the routine "setDirectKeys") is a subset of "tracers(:)"
92!     is extracted: the indexes are no longer valid for a subset (examples: iqParent, iqDescen).
93!     * If you need to convert a %key/%val pair into a direct-access key, add the corresponding line in "setDirectKeys".
94!==============================================================================================================================
95LOGICAL FUNCTION readTracersFiles(type_trac, fType, tracs) RESULT(lerr)
96!------------------------------------------------------------------------------------------------------------------------------
97  CHARACTER(LEN=*),             INTENT(IN)  :: type_trac              !--- List of components used
98  INTEGER,                      INTENT(OUT) :: fType                  !--- Type of input file found
99  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) :: tracs(:)
100  CHARACTER(LEN=maxlen),  ALLOCATABLE :: s(:), sections(:), trac_files(:)
101  CHARACTER(LEN=maxlen) :: str, fname, mesg
102  INTEGER               :: is, nsec, ierr, it, ntrac, ns, ip, ix
103  LOGICAL, ALLOCATABLE  :: ll(:), lGen3(:)
104!------------------------------------------------------------------------------------------------------------------------------
105  lerr = .FALSE.
106  modname = 'readTracersFiles'
107  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
108
109  !--- Required sections + corresponding files names (new style single section case)
110  IF(test(strParse(type_trac, '|', sections), lerr)) RETURN          !--- Parse "type_trac" list
111
112  nsec = SIZE(sections, DIM=1)
113  ALLOCATE(trac_files(nsec)); DO is=1, nsec; trac_files(is) = 'tracer_'//TRIM(sections(is))//'.def'; END DO
114
115  !--- LOOK AT AVAILABLE FILES
116  ll = .NOT.testFile(trac_files)
117  fType = 0
118  IF(.NOT.testFile('traceur.def') .AND. nsec==1) fType = 1           !--- OLD STYLE FILE
119  IF(.NOT.testFile('tracer.def'))                fType = 2           !--- NEW STYLE ; SINGLE  FILE, SEVERAL SECTIONS
120  IF(ALL(ll))                                    fType = 3           !--- NEW STYLE ; SEVERAL FILES, SINGLE SECTION USED
121  IF(ANY(ll) .AND. fType/=3) THEN                                    !--- MISSING FILES
122    IF(test(checkList(trac_files, .NOT.ll, 'Failed reading tracers description', 'files', 'missing'), lerr)) RETURN
123  END IF
124
125  !--- CHECK WHETHER type_trac AND FILE TYPE ARE COMPATIBLE
126  IF(test(fmsg('No multiple sections for the old format "traceur.def"', ll = SIZE(sections)>1 .AND. fType==1), lerr)) RETURN
127
128  !--- TELLS WHAT WAS IS ABOUT TO BE USED
129  IF (fmsg('No adequate tracers description file(s) found ; default values will be used',          modname, fType==0)) RETURN
130  CALL msg('Trying to read old-style tracers description file "traceur.def"',                      modname, fType==1)
131  CALL msg('Trying to read the new style multi-sections tracers description file "tracer.def"',    modname, fType==2)
132  CALL msg('Trying to read the new style single section tracers description files "tracer_*.def"', modname, fType==3)
133
134  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
135  SELECT CASE(fType)                         !--- Set %name, %genOName, %parent, %type, %phase, %component, %iGeneration, %keys
136  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
137    CASE(1)                                                               !=== OLD FORMAT "traceur.def"
138    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139      !--- OPEN THE "traceur.def" FILE
140      OPEN(90, FILE="traceur.def", FORM='formatted', STATUS='old', IOSTAT=ierr)
141
142      !--- GET THE TRACERS NUMBER
143      READ(90,'(i3)',IOSTAT=ierr)ntrac                               !--- Number of lines/tracers
144      IF(test(fmsg('Invalid format for "'//TRIM(fname)//'"', modname, ierr /= 0), lerr)) RETURN
145
146      !--- READ THE REMAINING LINES: <hadv> <vadv> <tracer> [<transporting fluid>]
147      ALLOCATE(tracs(ntrac))
148      DO it=1,ntrac                                                  !=== READ RAW DATA: loop on the line/tracer number
149        READ(90,'(a)',IOSTAT=ierr) str
150        IF(test(fmsg('Invalid format for "' //TRIM(fname)//'"', modname, ierr>0), lerr)) RETURN
151        IF(test(fmsg('Not enough lines in "'//TRIM(fname)//'"', modname, ierr<0), lerr)) RETURN
152        ll = strParse(str, ' ', s, n=ns)
153        CALL msg('This file is for air tracers only',           modname, ns == 3 .AND. it == 1)
154        CALL msg('This files specifies the transporting fluid', modname, ns == 4 .AND. it == 1)
155        tracs(it)%name   = old2newName(s(3), ip)                     !--- Set %name:   name   of the tracer
156        tracs(it)%parent = tran0                                     !--- Default transporting fluid name
157        IF(ns == 4) tracs(it)%parent = old2newName(s(4))             !--- Set %parent: parent of the tracer
158        tracs(it)%phase = known_phases(ip:ip)                        !--- Set %phase:  tracer phase (default: "g"azeous)
159        tracs(it)%component = TRIM(type_trac)                        !--- Set %component: model component name
160        tracs(it)%keys%key = ['hadv', 'vadv']                        !--- Set %keys%key
161        tracs(it)%keys%val = s(1:2)                                  !--- Set %keys%val
162      END DO
163      CLOSE(90)
164      CALL setGeneration(tracs)                                      !--- Set %iGeneration and %gen0Name
165      WHERE(tracs%iGeneration == 2) tracs%type = 'tag'               !--- Set %type:        'tracer' or 'tag'
166      IF(test(checkTracers(tracs, fname, fname), lerr)) RETURN       !--- Detect orphans and check phases
167      IF(test(checkUnique (tracs, fname, fname), lerr)) RETURN       !--- Detect repeated tracers
168      CALL sortTracers  (tracs)                                      !--- Sort the tracers
169      tracs(:)%keys%name = tracs(:)%name                             !--- Copy tracers names in keys components
170    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
171    CASE(2); IF(test(feedDBase(["tracer.def"], [type_trac], modname), lerr)) RETURN !=== SINGLE   FILE, MULTIPLE SECTIONS
172    !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
173    CASE(3); IF(test(feedDBase(  trac_files  ,  sections,   modname), lerr)) RETURN !=== MULTIPLE FILES, SINGLE  SECTION
174  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
175  END SELECT
176  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
177
178
179  IF(ALL([2,3] /= fType)) RETURN
180
181  IF(nsec  == 1) THEN;
182    tracs = dBase(1)%trac
183  ELSE IF(tracs_merge) THEN
184    CALL msg('The multiple required sections will be MERGED.',    modname)
185    IF(test(mergeTracers(dBase, tracs), lerr)) RETURN
186  ELSE
187    CALL msg('The multiple required sections will be CUMULATED.', modname)
188    IF(test(cumulTracers(dBase, tracs), lerr)) RETURN
189  END IF
190  WHERE(tracs%gen0Name(1:3) /= 'H2O') tracs%isInPhysics=.TRUE.       !--- Set %isInPhysics: passed to physics
191  CALL setDirectKeys(tracs)                                          !--- Set %iqParent, %iqDescen, %nqDescen, %nqChilds
192END FUNCTION readTracersFiles
193!==============================================================================================================================
194
195!==============================================================================================================================
196LOGICAL FUNCTION feedDBase(fnames, snames, modname) RESULT(lerr)
197! Purpose: Read the sections "snames(is)" (pipe-separated list) from each "fnames(is)"
198!   file and create the corresponding tracers set descriptors in the database "dBase":
199! * dBase(id)%name                : section name
200! * dBase(id)%trac(:)%name        : tracers names
201! * dBase(id)%trac(it)%keys%key(:): names  of keys associated to tracer dBase(id)%trac(it)%name
202! * dBase(id)%trac(it)%keys%val(:): values of keys associated to tracer dBase(id)%trac(it)%name
203!------------------------------------------------------------------------------------------------------------------------------
204  CHARACTER(LEN=*), INTENT(IN)  :: fnames(:)                         !--- Files names
205  CHARACTER(LEN=*), INTENT(IN)  :: snames(:)                         !--- Pipe-deparated list of sections (one list each file)
206  CHARACTER(LEN=*), INTENT(IN)  :: modname                           !--- Calling routine name
207  INTEGER,  ALLOCATABLE :: ndb(:)                                    !--- Number of sections for each file
208  INTEGER,  ALLOCATABLE :: ixf(:)                                    !--- File index for each section of the expanded list
209  LOGICAL,  ALLOCATABLE :: lTg(:)                                    !--- Tagging tracers mask
210  CHARACTER(LEN=maxlen) :: fnm, snm
211  INTEGER               :: idb, i
212  LOGICAL :: ll
213!------------------------------------------------------------------------------------------------------------------------------
214  !=== READ THE REQUIRED SECTIONS
215  ll = strCount(snames, '|', ndb)                                    !--- Number of sections for each file
216  ALLOCATE(ixf(SUM(ndb)))
217  DO i=1, SIZE(fnames)                                               !--- Set %name, %keys
218    IF(test(readSections(fnames(i), snames(i), 'default'), lerr)) RETURN
219    ixf(1+SUM(ndb(1:i-1)):SUM(ndb(1:i))) = i                         !--- File index for each section of the expanded list
220  END DO
221  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
222  DO idb=1,SIZE(dBase)                                               !--- LOOP ON THE LOADED SECTIONS
223  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
224    fnm = fnames(ixf(idb)); snm = dBase(idb)%name                    !--- FILE AND SECTION NAMES
225    lerr = ANY([(dispTraSection('RAW CONTENT OF SECTION "'//TRIM(snm)//'"', snm, modname), idb=1, SIZE(dBase))])
226    IF(test(expandSection(dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- EXPAND NAMES ;  set %parent, %type, %component
227    CALL setGeneration   (dBase(idb)%trac)                           !---                 set %iGeneration,   %genOName
228    IF(test(checkTracers (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK ORPHANS AND PHASES
229    IF(test(checkUnique  (dBase(idb)%trac, snm, fnm), lerr)) RETURN  !--- CHECK TRACERS UNIQUENESS
230    CALL expandPhases    (dBase(idb)%trac)                           !--- EXPAND PHASES ; set %phase
231    CALL sortTracers     (dBase(idb)%trac)                           !--- SORT TRACERS
232    lerr = ANY([(dispTraSection('EXPANDED CONTENT OF SECTION "'//TRIM(snm)//'"', snm, modname), idb=1, SIZE(dBase))])
233  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
234  END DO
235  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
236END FUNCTION feedDBase
237!------------------------------------------------------------------------------------------------------------------------------
238
239!------------------------------------------------------------------------------------------------------------------------------
240LOGICAL FUNCTION readSections(fnam,snam,defName) RESULT(lerr)
241!------------------------------------------------------------------------------------------------------------------------------
242  CHARACTER(LEN=*),           INTENT(IN) :: fnam                     !--- File name
243  CHARACTER(LEN=*),           INTENT(IN) :: snam                     !--- Pipe-separated sections list
244  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: defName                  !--- Special section (default values) name
245!------------------------------------------------------------------------------------------------------------------------------
246  TYPE(dataBase_type),   ALLOCATABLE :: tdb(:)
247  CHARACTER(LEN=maxlen), ALLOCATABLE :: sec(:)
248  INTEGER,               ALLOCATABLE ::  ix(:)
249  INTEGER :: n0, idb, ndb, i, j
250  LOGICAL :: ll
251!------------------------------------------------------------------------------------------------------------------------------
252  n0 = SIZE(dBase) + 1                                               !--- Index for next entry in the database
253  CALL readSections_all()                                            !--- Read all the sections of file "fnam"
254  ndb= SIZE(dBase)                                                   !--- Current number of sections in the database
255  IF(PRESENT(defName)) THEN                                          !--- Add default values to all the tracers
256    DO idb=n0,ndb; CALL addDefault(dBase(idb)%trac, defName); END DO !--- and remove the virtual tracer "defName"
257  END IF
258  ll = strParse(snam, '|', keys = sec)                               !--- Requested sections names
259  ix = strIdx(dBase(:)%name, sec(:))                                 !--- Indexes of requested sections in database
260  IF(test(checkList(sec, ix == 0, 'In file "'//TRIM(fnam)//'"','section(s)', 'missing'), lerr)) RETURN
261  tdb = dBase(:); dBase = [tdb(1:n0-1),tdb(PACK(ix, MASK=ix/=0))]    !--- Keep requested sections only
262
263CONTAINS
264
265!------------------------------------------------------------------------------------------------------------------------------
266SUBROUTINE readSections_all()
267!------------------------------------------------------------------------------------------------------------------------------
268  CHARACTER(LEN=maxlen), ALLOCATABLE ::  s(:), v(:)
269  TYPE(trac_type),       ALLOCATABLE :: tt(:)
270  TYPE(trac_type)       :: tmp
271  CHARACTER(LEN=1024)   :: str, str2
272  CHARACTER(LEN=maxlen) :: secn
273  INTEGER               :: ierr, n
274!------------------------------------------------------------------------------------------------------------------------------
275  IF(.NOT.ALLOCATED(dBase)) ALLOCATE(dBase(0))
276  OPEN(90, FILE=fnam, FORM='formatted', STATUS='old')
277  DO; str=''
278    DO
279      READ(90,'(a)', IOSTAT=ierr)str2                                !--- Read a full line
280      str=TRIM(str)//' '//TRIM(str2)                                 !--- Append "str" with the current line
281      n=LEN_TRIM(str); IF(n == 0) EXIT                               !--- Empty line (probably end of file)
282      IF(IACHAR(str(n:n))  /= 92) EXIT                               !--- No "\" continuing line symbol found => end of line
283      str = str(1:n-1)                                               !--- Remove the "\" continuing line symbol
284    END DO
285    str = ADJUSTL(str)                                               !--- Remove the front space
286    IF(ierr    /= 0 ) EXIT                                           !--- Finished: error or end of file
287    IF(str(1:1)=='#') CYCLE                                          !--- Skip comments lines
288    CALL removeComment(str)                                          !--- Skip comments at the end of a line
289    IF(str     == '') CYCLE                                          !--- Skip empty line (probably at the end of the file)
290    IF(str(1:1)=='&') THEN                                           !=== SECTION HEADER LINE
291      ndb  = SIZE(dBase)                                             !--- Number of sections so far
292      secn = str(2:LEN_TRIM(str))//' '                               !--- Current section name
293      IF(ANY(dBase(:)%name == secn)) CYCLE                           !--- Already known section
294      IF(secn(1:7) == 'version')     CYCLE                           !--- Skip the "version" special section
295      ndb = ndb + 1                                                  !--- Extend database
296      ALLOCATE(tdb(ndb))
297      tdb(1:ndb-1)  = dBase
298      tdb(ndb)%name = secn
299      ALLOCATE(tdb(ndb)%trac(0))
300      CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
301    ELSE                                                             !=== TRACER LINE
302      ll = strParse(str,' ', keys = s, vals = v, n = n)              !--- Parse <key>=<val> pairs
303      tt = dBase(ndb)%trac(:)
304      tmp%name = s(1); tmp%keys = keys_type(s(1), s(2:n), v(2:n))    !--- Set %name and %keys
305      dBase(ndb)%trac = [tt(:), tmp]
306      DEALLOCATE(tt)
307!      dBase(ndb)%trac = [dBase(ndb)%trac(:), tra(name=s(1), keys=keys_type(s(1), s(2:n), v(2:n)))]
308    END IF
309  END DO
310  CLOSE(90)
311
312END SUBROUTINE readSections_all
313!------------------------------------------------------------------------------------------------------------------------------
314
315END FUNCTION readSections
316!==============================================================================================================================
317
318
319!==============================================================================================================================
320SUBROUTINE addDefault(t, defName)
321!------------------------------------------------------------------------------------------------------------------------------
322! Purpose: Add the keys from virtual tracer named "defName" (if any) and remove this virtual tracer.
323!------------------------------------------------------------------------------------------------------------------------------
324  TYPE(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
325  CHARACTER(LEN=*),                     INTENT(IN)    :: defName
326  INTEGER :: jd, it, k
327  TYPE(keys_type), POINTER :: ky
328  TYPE(trac_type), ALLOCATABLE :: tt(:)
329  jd = strIdx(t(:)%name, defName)
330  IF(jd == 0) RETURN
331  ky => t(jd)%keys
332  DO k = 1, SIZE(ky%key)                                             !--- Loop on the keys of the tracer named "defName"
333!   CALL addKey_m(ky%key(k), ky%val(k), t(:)%keys)                   !--- Add key to all the tracers (no overwriting)
334    DO it = 1, SIZE(t); CALL addKey_1(ky%key(k), ky%val(k), t(it)%keys); END DO
335  END DO
336  tt = [t(1:jd-1),t(jd+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
337END SUBROUTINE addDefault
338!==============================================================================================================================
339
340!==============================================================================================================================
341SUBROUTINE subDefault(t, defName, lSubLocal)
342!------------------------------------------------------------------------------------------------------------------------------
343! Purpose: Substitute the keys from virtual tracer named "defName" (if any) and remove this virtual tracer.
344!          Substitute the keys locally (for the current tracer) if the flag "lSubLocal" is .TRUE.
345!------------------------------------------------------------------------------------------------------------------------------
346  TYPE(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: t(:)
347  CHARACTER(LEN=*),                     INTENT(IN)    :: defName
348  LOGICAL,                              INTENT(IN)    :: lSubLocal
349  INTEGER :: i0, it, ik
350  TYPE(keys_type), POINTER     :: k0, ky
351  TYPE(trac_type), ALLOCATABLE :: tt(:)
352  i0 = strIdx(t(:)%name, defName)
353  IF(i0 == 0) RETURN
354  k0 => t(i0)%keys
355  DO it = 1, SIZE(t); IF(it == i0) CYCLE                             !--- Loop on the tracers
356    ky => t(it)%keys
357
358    !--- Substitute in the values of <key>=<val> pairs the keys defined in the virtual tracer "defName"
359    DO ik = 1, SIZE(k0%key); CALL strReplace(ky%val, k0%key(ik), k0%val(ik), .TRUE.); END DO
360
361    IF(.NOT.lSubLocal) CYCLE
362    !--- Substitute in the values of <key>=<val> pairs the keys defined locally (in the current tracer)
363    DO ik = 1, SIZE(ky%key); CALL strReplace(ky%val, ky%key(ik), ky%val(ik), .TRUE.); END DO
364  END DO
365  tt = [t(1:i0-1),t(i0+1:SIZE(t))]; CALL MOVE_ALLOC(FROM=tt, TO=t)   !--- Remove the virtual tracer named "defName"
366
367END SUBROUTINE subDefault
368!==============================================================================================================================
369
370
371!==============================================================================================================================
372LOGICAL FUNCTION expandSection(tr, sname, fname) RESULT(lerr)
373!------------------------------------------------------------------------------------------------------------------------------
374! Purpose: Expand tracers and parents lists in the tracers descriptor "tra".
375! Note:  * The following keys are expanded, so are accessible explicitely using "%" operator: "parent" "type".
376!        * Default values are provided for these keys because they are necessary.
377!------------------------------------------------------------------------------------------------------------------------------
378  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)                 !--- Tracer derived type vector
379  CHARACTER(LEN=*),             INTENT(IN)    :: sname
380  CHARACTER(LEN=*), OPTIONAL,   INTENT(IN)    :: fname
381  TYPE(trac_type),       ALLOCATABLE :: ttr(:)
382  CHARACTER(LEN=maxlen), ALLOCATABLE :: ta(:), pa(:)
383  CHARACTER(LEN=maxlen) :: msg1, modname
384  INTEGER :: it, nt, iq, nq, itr, ntr, ipr, npr, i
385  LOGICAL :: ll
386  modname = 'expandSection'
387  lerr = .FALSE.
388  nt = SIZE(tr)
389  nq = 0
390  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
391  DO it = 1, nt    !=== GET TRACERS NB AFTER EXPANSION + NEEDED KEYS (name, parent, type)
392  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
393    !--- Extract useful keys: parent name, type, component name
394    tr(it)%parent    = fgetKey(it, 'parent', tr(:)%keys,  tran0  )
395    tr(it)%type      = fgetKey(it, 'type'  , tr(:)%keys, 'tracer')
396    tr(it)%component = sname
397
398    !--- Determine the number of tracers and parents ; coherence checking
399    ll = strCount(tr(it)%name,   ',', ntr)
400    ll = strCount(tr(it)%parent, ',', npr)
401
402    !--- Tagging tracers only can have multiple parents
403    IF(test(npr/=1 .AND. TRIM(tr(it)%type)/='tag', lerr)) THEN
404      msg1 = 'Check section "'//TRIM(sname)//'"'
405      IF(PRESENT(fname)) msg1=TRIM(msg1)//' in file "'//TRIM(fname)//'"'
406      CALL msg(TRIM(msg1)//': "'//TRIM(tr(it)%name)//'" has several parents but is not a tag', modname); RETURN
407    END IF
408    nq = nq + ntr*npr                 
409  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
410  END DO
411  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
412  CALL delKey(['parent','type  '], tr)
413
414  ALLOCATE(ttr(nq))
415  iq = 1
416  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
417  DO it = 1, nt                                                      !=== EXPAND TRACERS AND PARENTS NAMES LISTS
418  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
419    ll = strParse(tr(it)%name,   ',', ta, n=ntr)                     !--- Number of tracers
420    ll = strParse(tr(it)%parent, ',', pa, n=npr)                     !--- Number of parents
421    DO ipr=1,npr                                                     !--- Loop on parents list elts
422      DO itr=1,ntr                                                   !--- Loop on tracers list elts
423        i = iq+itr-1+(ipr-1)*ntr
424        ttr(i)%name   = TRIM(ta(itr))
425        ttr(i)%parent = TRIM(pa(ipr))
426        ttr(i)%keys%name = ta(itr)
427        ttr(i)%keys%key  = tr(it)%keys%key
428        ttr(i)%keys%val  = tr(it)%keys%val
429!        ttr(i)%keys   = keys_type(ta(itr), tr(it)%keys%key, tr(it)%keys%val)
430      END DO
431    END DO
432    ttr(iq:iq+ntr*npr-1)%type      = tr(it)%type                     !--- Duplicating type
433    ttr(iq:iq+ntr*npr-1)%component = tr(it)%component                !--- Duplicating type
434    iq = iq + ntr*npr
435  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
436  END DO
437  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
438  DEALLOCATE(ta,pa)
439  CALL MOVE_ALLOC(FROM=ttr, TO=tr)
440
441END FUNCTION expandSection
442!==============================================================================================================================
443
444!==============================================================================================================================
445SUBROUTINE setGeneration(tr)
446!------------------------------------------------------------------------------------------------------------------------------
447! Purpose: Determine, for each tracer of "tr(:)":
448!   * %iGeneration: the generation number
449!   * %gen0Name:    the generation 0 ancestor name
450!------------------------------------------------------------------------------------------------------------------------------
451  TYPE(trac_type),     INTENT(INOUT) :: tr(:)                        !--- Tracer derived type vector
452  INTEGER                            :: iq, nq, ig
453  LOGICAL,               ALLOCATABLE :: lg(:)
454  CHARACTER(LEN=maxlen), ALLOCATABLE :: prn(:)
455!------------------------------------------------------------------------------------------------------------------------------
456  tr(:)%iGeneration = -1                                             !--- error if -1
457  nq = SIZE(tr, DIM=1)                                               !--- Number of tracers lines
458  lg = tr(:)%parent == tran0                                         !--- Flag for generation 0 tracers
459  WHERE(lg) tr(:)%iGeneration = 0                                    !--- Generation 0 tracers
460
461  !=== Determine generation for each tracer
462  ig=-1; prn = [tran0]
463  DO                                                                 !--- Update current generation flag
464    IF(ig/=-1) prn = PACK( tr(:)%name, MASK=tr(:)%iGeneration == ig)
465    lg(:) = [(ANY(prn(:) == tr(iq)%parent), iq=1, nq)]               !--- Current generation tracers flag
466    IF( ALL( .NOT. lg ) ) EXIT                                       !--- Empty current generation
467    ig = ig+1; WHERE(lg) tr(:)%iGeneration = ig
468  END DO
469  tr(:)%gen0Name = ancestor(tr)                                      !--- First generation ancestor name
470
471END SUBROUTINE setGeneration
472!==============================================================================================================================
473
474!==============================================================================================================================
475LOGICAL FUNCTION checkTracers(tr, sname, fname) RESULT(lerr)
476!------------------------------------------------------------------------------------------------------------------------------
477! Purpose:
478!   * check for orphan tracers (without known parent)
479!   * check wether the phases are known or not ("g"aseous, "l"iquid or "s"olid so far)
480!------------------------------------------------------------------------------------------------------------------------------
481  TYPE(trac_type),            INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
482  CHARACTER(LEN=*),           INTENT(IN) :: sname                    !--- Section name
483  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
484  CHARACTER(LEN=maxlen) :: mesg
485  CHARACTER(LEN=maxlen) :: bp(SIZE(tr, DIM=1)), pha                  !--- Bad phases list, phases of current tracer
486  CHARACTER(LEN=1) :: p
487  INTEGER :: ip, np, iq, nq
488!------------------------------------------------------------------------------------------------------------------------------
489  nq = SIZE(tr,DIM=1)                                                !--- Number of tracers lines
490  mesg = 'Check section "'//TRIM(sname)//'"'
491  IF(PRESENT(fname)) mesg=TRIM(mesg)//' in file "'//TRIM(fname)//'"'
492
493  !=== CHECK FOR ORPHAN TRACERS
494  IF(test(checkList(tr%name, tr%iGeneration==-1, mesg, 'tracers', 'orphan'), lerr)) RETURN
495
496  !=== CHECK PHASES
497  DO iq=1,nq; IF(tr(iq)%iGeneration/=0) CYCLE                        !--- Generation O only is checked
498    pha = fgetKey(iq, 'phases', tr(:)%keys, 'g')                     !--- Phases
499    np = LEN_TRIM(pha); bp(iq)=' '
500    DO ip=1,np; p = pha(ip:ip); IF(INDEX(known_phases,p)==0) bp(iq) = TRIM(bp(iq))//p; END DO
501    IF(TRIM(bp(iq)) /= '') bp(iq) = TRIM(tr(iq)%name)//': '//TRIM(bp(iq))
502  END DO
503  lerr = checkList(bp, tr%iGeneration==0 .AND. bp/='', mesg, 'tracers phases', 'unknown')
504END FUNCTION checkTracers
505!==============================================================================================================================
506
507!==============================================================================================================================
508LOGICAL FUNCTION checkUnique(tr, sname, fname) RESULT(lerr)
509!------------------------------------------------------------------------------------------------------------------------------
510! Purpose: Make sure that tracers are not repeated.
511!------------------------------------------------------------------------------------------------------------------------------
512  TYPE(trac_type),            INTENT(IN) :: tr(:)                    !--- Tracer derived type vector
513  CHARACTER(LEN=*),           INTENT(IN) :: sname                    !--- Section name
514  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fname                    !--- File name
515!------------------------------------------------------------------------------------------------------------------------------
516  INTEGER :: ip, np, iq, nq, k
517  LOGICAL, ALLOCATABLE  :: ll(:)
518  CHARACTER(LEN=maxlen) :: mesg, tnam, tdup(SIZE(tr,DIM=1))
519  CHARACTER(LEN=1)      :: p
520!------------------------------------------------------------------------------------------------------------------------------
521  mesg = 'Check section "'//TRIM(sname)//'"'
522  IF(PRESENT(fname)) mesg=TRIM(mesg)//' in file "'//TRIM(fname)//'"'
523  nq=SIZE(tr,DIM=1); lerr=.FALSE.                                    !--- Number of lines ; error flag
524  tdup(:) = ''
525  DO iq=1,nq; IF(tr(iq)%type == 'tag') CYCLE                         !--- Tags can be repeated
526    tnam = TRIM(tr(iq)%name)
527    ll = tr(:)%name==TRIM(tnam)                                      !--- Mask for current tracer name
528    IF(COUNT(ll)==1 ) CYCLE                                          !--- Tracer is not repeated
529    IF(tr(iq)%iGeneration>0) THEN
530      tdup(iq) = tnam                                                !--- gen>0: MUST be unique
531    ELSE
532      DO ip=1,nphases; p=known_phases(ip:ip)                         !--- Loop on known phases
533        !--- Number of appearances of the current tracer with known phase "p"
534        np = COUNT( PACK( [(INDEX(fgetKey(k, 'phases', tr(:)%keys, 'g'),p), k=1, nq)] /=0 , MASK=ll ) )
535        IF(np <=1) CYCLE
536        tdup(iq) = TRIM(tdup(iq))//TRIM(phases_names(ip))
537        IF(ANY(tdup(1:iq-1) == tdup(iq))) tdup(iq)=''                !--- Avoid repeating same messages
538      END DO
539    END IF
540    IF(tdup(iq) /= '') tdup(iq)=TRIM(tnam)//' in '//TRIM(tdup(iq))//' phase(s)'
541  END DO
542  lerr = checkList(tdup, tdup/='', mesg, 'tracers', 'duplicated')
543END FUNCTION checkUnique
544!==============================================================================================================================
545
546!==============================================================================================================================
547SUBROUTINE expandPhases(tr)
548!------------------------------------------------------------------------------------------------------------------------------
549! Purpose: Expand the phases in the tracers descriptor "tr". Phases are not repeated for a tracer, thanks to "checkUnique".
550!------------------------------------------------------------------------------------------------------------------------------
551  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)               !--- Tracer derived type vector
552!------------------------------------------------------------------------------------------------------------------------------
553  TYPE(trac_type), ALLOCATABLE :: ttr(:)
554  INTEGER,   ALLOCATABLE ::  i0(:)
555  CHARACTER(LEN=maxlen)  :: nam, pha, trn
556  CHARACTER(LEN=1) :: p
557  INTEGER :: ip, np, iq, jq, nq, it, nt, nc, i, n
558  LOGICAL :: lTg, lEx
559!------------------------------------------------------------------------------------------------------------------------------
560  nq = SIZE(tr, DIM=1)
561  nt = 0
562  DO iq = 1, nq                                                      !--- GET THE NUMBER OF TRACERS
563    IF(tr(iq)%iGeneration /= 0) CYCLE                                !--- Only deal with generation 0 tracers
564    nc = COUNT(tr(:)%gen0Name==tr(iq)%name .AND. tr%iGeneration/=0)  !--- Number of childs of tr(iq)
565    tr(iq)%phase = fgetKey(iq, 'phases', tr(:)%keys)                 !--- Phases list      of tr(iq)
566    np = LEN_TRIM(tr(iq)%phase)                                      !--- Number of phases of tr(iq)
567    nt = nt + (1+nc) * np                                            !--- Number of tracers after expansion
568  END DO
569  ALLOCATE(ttr(nt))                                                  !--- Version  of "tr" after phases expansion
570  it = 1                                                             !--- Current "ttr(:)" index
571  DO iq = 1, nq                                                      !--- Loop on "tr(:)" indexes
572    lTg = tr(iq)%type=='tag'                                         !--- Current tracer is a tag
573    i0 = strFind(tr(:)%name, TRIM(tr(iq)%gen0Name), n)               !--- Indexes of first generation ancestor copies
574    np = SUM([( LEN_TRIM(tr(i0(i))%phase),i=1,n )], 1)               !--- Number of phases for current tracer tr(iq)
575    lEx = np>1                                                       !--- Phase suffix only required if phases number is > 1
576    IF(lTg) lEx = lEx .AND. tr(iq)%iGeneration>0                     !--- No phase suffix for generation 0 tags
577    DO i=1,n                                                         !=== LOOP ON GENERATION 0 ANCESTORS
578      jq = i0(i)                                                     !--- tr(jq): ith tracer with same gen 0 ancestor as tr(iq)
579      IF(tr(iq)%iGeneration==0) jq=iq                                !--- Generation 0: count the current tracer phases only
580      pha = tr(jq)%phase                                             !--- Phases list for tr(jq)
581      DO ip = 1, LEN_TRIM(pha)                                       !=== LOOP ON PHASES LISTS
582        p = pha(ip:ip)
583        trn = TRIM(tr(iq)%name); nam = trn                           !--- Tracer name (regular case)
584        IF(lTg) nam = TRIM(tr(iq)%parent)                            !--- Parent name (tagging case)
585        IF(lEx) nam = addPhase(nam, p )                              !--- Phase extension needed
586        IF(lTg) nam = TRIM(nam)//'_'//TRIM(trn)                      !--- <parent>_<name> for tags
587        ttr(it) = tr(iq)                                             !--- Same <key>=<val> pairs
588        ttr(it)%name      = TRIM(nam)                                !--- Name with possibly phase suffix
589        ttr(it)%keys%name = TRIM(nam)                                !--- Name inside the keys decriptor
590        ttr(it)%phase     = p                                        !--- Single phase entry
591        IF(lEx .AND. tr(iq)%iGeneration>0) THEN
592          ttr(it)%parent   = addPhase(ttr(it)%parent,   p)
593          ttr(it)%gen0Name = addPhase(ttr(it)%gen0Name, p)
594        END IF
595        it = it+1
596      END DO
597      IF(tr(iq)%iGeneration==0) EXIT                                 !--- Break phase loop for gen 0
598    END DO
599  END DO
600  CALL MOVE_ALLOC(FROM=ttr, TO=tr)
601  CALL delKey(['phases'],tr)                                         !--- Remove few keys entries
602
603END SUBROUTINE expandPhases
604!==============================================================================================================================
605
606!==============================================================================================================================
607SUBROUTINE sortTracers(tr)
608!------------------------------------------------------------------------------------------------------------------------------
609! Purpose: Sort tracers:
610!  * Put water at the beginning of the vector, in the "known_phases" order.
611!  * lGrowGen == T: in ascending generations numbers.
612!  * lGrowGen == F: tracer + its childs sorted by growing generation, one after the other.
613!   TO BE ADDED IF NECESSARY: HIGHER MOMENTS AT THE END
614!------------------------------------------------------------------------------------------------------------------------------
615  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: tr(:)         !--- Tracer derived type vector
616!------------------------------------------------------------------------------------------------------------------------------
617  TYPE(trac_type), ALLOCATABLE        :: tr2(:)
618  INTEGER,         ALLOCATABLE        :: iy(:), iz(:)
619  INTEGER                             :: ig, ng, iq, jq, ip, nq, n, ix(SIZE(tr)), k
620  INTEGER                             :: it
621!  tr2 is introduced in order to cope with a bug in gfortran 4.8.5 compiler
622!------------------------------------------------------------------------------------------------------------------------------
623  nq = SIZE(tr)
624  DO ip = nphases, 1, -1
625    iq = strIdx(tr(:)%name, addPhase('H2O', ip))
626    IF(iq == 0) CYCLE
627    tr2 = tr(:)
628    tr = [tr2(iq), tr2(1:iq-1), tr2(iq+1:nq)]
629  END DO
630  IF(lSortByGen) THEN
631    iq = 1
632    ng = MAXVAL(tr(:)%iGeneration, MASK=.TRUE., DIM=1)               !--- Number of generations
633    DO ig = 0, ng                                                    !--- Loop on generations
634      iy = PACK([(k, k=1, nq)], MASK=tr(:)%iGeneration==ig)          !--- Generation ig tracers indexes
635      n = SIZE(iy)
636      ix(iq:iq+n-1) = iy                                             !--- Stack growing generations idxs
637      iq = iq + n
638    END DO
639  ELSE
640    iq = 1
641    DO jq = 1, nq                                                    !--- Loop on generation 0 tracers
642      IF(tr(jq)%iGeneration /= 0) CYCLE                              !--- Skip generations /= 0
643      ix(iq) = jq                                                    !--- Generation 0 ancestor index first
644      iq = iq + 1                                                    !--- Next "iq" for next generations tracers
645      iy = strFind(tr(:)%gen0Name, TRIM(tr(jq)%name))                !--- Indexes of "tr(jq)" childs in "tr(:)"
646      ng = MAXVAL(tr(iy)%iGeneration, MASK=.TRUE., DIM=1)            !--- Number of generations of the "tr(jq)" family
647      DO ig = 1, ng                                                  !--- Loop   on generations of the "tr(jq)" family
648        iz = find(tr(iy)%iGeneration, ig, n)                         !--- Indexes of the tracers "tr(iy(:))" of generation "ig"
649        ix(iq:iq+n-1) = iy(iz)                                       !--- Same indexes in "tr(:)"
650        iq = iq + n
651      END DO
652    END DO
653  END IF
654  tr = tr(ix)                                                        !--- Reorder the tracers
655END SUBROUTINE sortTracers
656!==============================================================================================================================
657
658!==============================================================================================================================
659LOGICAL FUNCTION mergeTracers(sections, tr) RESULT(lerr)
660  TYPE(dataBase_type),  TARGET, INTENT(IN)  :: sections(:)
661  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) ::       tr(:)
662  TYPE(trac_type), POINTER ::   t1(:),   t2(:)
663  INTEGER,     ALLOCATABLE :: ixct(:), ixck(:)
664  INTEGER :: is, k1, k2, nk2, i1, i2, nt2
665  CHARACTER(LEN=maxlen) :: s1, v1, v2, tnam, knam, modname
666  modname = 'mergeTracers'
667  lerr = .FALSE.
668  t1 => sections(1)%trac(:)                                          !--- Alias: first tracers section
669  tr = t1
670  !----------------------------------------------------------------------------------------------------------------------------
671  DO is=2,SIZE(sections)                                             !=== SEVERAL SECTIONS: MERGE THEM
672  !----------------------------------------------------------------------------------------------------------------------------
673    t2  => sections(is)%trac(:)                                      !--- Alias: current tracers section
674    nt2  = SIZE(t2(:), DIM=1)                                        !--- Number of tracers in section
675    ixct = strIdx(t1(:)%name, t2(:)%name)                            !--- Indexes of common tracers
676    tr = [tr, PACK(t2, MASK= ixct==0)]                               !--- Append with new tracers
677    IF( ALL(ixct == 0) ) CYCLE                                       !--- No common tracers => done
678    CALL msg('Tracers defined in previous sections and duplicated in "'//TRIM(sections(is)%name)//'":', modname)
679    CALL msg(t1(PACK(ixct, MASK = ixct/=0))%name, modname, nmax=128) !--- Display duplicates (the 128 first at most)
680    !--------------------------------------------------------------------------------------------------------------------------
681    DO i2=1,nt2; tnam = TRIM(t2(i2)%name)                            !=== LOOP ON COMMON TRACERS
682    !--------------------------------------------------------------------------------------------------------------------------
683      i1 = ixct(i2); IF(i1 == 0) CYCLE                               !--- Idx in t1(:) ; skip new tracers
684
685      !=== CHECK WETHER ESSENTIAL KEYS ARE IDENTICAL OR NOT
686      s1=' of "'//TRIM(tnam)//'" in "'//TRIM(sections(is)%name)//'" not matching previous value'
687     
688      IF(test(fmsg('Parent name'//TRIM(s1), modname, t1(i1)%parent      /= t2(i2)%parent),      lerr)) RETURN
689      IF(test(fmsg('Type'       //TRIM(s1), modname, t1(i1)%type        /= t2(i2)%type),        lerr)) RETURN
690      IF(test(fmsg('Generation' //TRIM(s1), modname, t1(i1)%iGeneration /= t2(i2)%iGeneration), lerr)) RETURN
691
692      !=== APPEND <key>=<val> PAIRS NOT PREVIOULSLY DEFINED
693      nk2  = SIZE(t2(i2)%keys%key(:))                                !--- Keys number in current section
694      ixck = strIdx(t1(i1)%keys%key(:), t2(i2)%keys%key(:))          !--- Common keys indexes
695
696      !=== APPEND NEW KEYS
697      tr(i1)%keys%key = [ tr(i1)%keys%key, PACK(tr(i2)%keys%key, MASK = ixck==0)]
698      tr(i1)%keys%val = [ tr(i1)%keys%val, PACK(tr(i2)%keys%val, MASK = ixck==0)]
699
700      !--- KEEP TRACK OF THE COMPONENTS NAMES
701      tr(i1)%component = TRIM(tr(i1)%component)//','//TRIM(tr(i2)%component)
702
703      !--- SELECT COMMON TRACERS WITH DIFFERING KEYS VALUES (PREVIOUS VALUE IS KEPT)
704      DO k2=1,nk2
705        k1 = ixck(k2); IF(k1 == 0) CYCLE
706        IF(t1(i1)%keys%val(k1) == t2(i2)%keys%val(k2)) ixck(k2)=0
707      END DO
708      IF(ALL(ixck==0)) CYCLE                                         !--- No identical keys with /=values
709
710      !--- DISPLAY INFORMATION: OLD VALUES ARE KEPT FOR THE KEYS FOUND IN PREVIOUS AND CURRENT SECTIONS
711      CALL msg('Key(s)'//TRIM(s1), modname)
712      DO k2 = 1, nk2                                                 !--- Loop on keys found in both t1(:) and t2(:)
713        knam = t2(i2)%keys%key(k2)                                   !--- Name of the current key
714        k1 = ixck(k2)                                                !--- Corresponding index in t1(:)
715        IF(k1 == 0) CYCLE                                            !--- New keys are skipped
716        v1 = t1(i1)%keys%val(k1); v2 = t2(i2)%keys%val(k2)           !--- Key values in t1(:) and t2(:)
717        CALL msg(' * '//TRIM(knam)//'='//TRIM(v2)//' ; previous value kept:'//TRIM(v1), modname)
718      END DO
719      !------------------------------------------------------------------------------------------------------------------------
720    END DO
721    !--------------------------------------------------------------------------------------------------------------------------
722  END DO
723  CALL sortTracers(tr)
724
725END FUNCTION mergeTracers
726!==============================================================================================================================
727
728!==============================================================================================================================
729LOGICAL FUNCTION cumulTracers(sections, tr) RESULT(lerr)
730  TYPE(dataBase_type),  TARGET, INTENT(IN)  :: sections(:)
731  TYPE(trac_type), ALLOCATABLE, INTENT(OUT) ::       tr(:)
732  TYPE(trac_type), POINTER     :: t1(:), t2(:)
733  INTEGER,   ALLOCATABLE :: nt(:)
734  CHARACTER(LEN=maxlen)  :: tnam, tnam_new
735  INTEGER :: iq, nq, is, ns, nsec
736  lerr = .FALSE.                                                     !--- Can't fail ; kept to match "mergeTracer" interface.
737  nsec =  SIZE(sections)
738  tr = [(      sections(is)%trac(:) , is=1, nsec )]                  !--- Concatenated tracers vector
739  nt = [( SIZE(sections(is)%trac(:)), is=1, nsec )]                  !--- Number of tracers in each section
740  !----------------------------------------------------------------------------------------------------------------------------
741  DO is=1, nsec                                                      !=== LOOP ON SECTIONS
742  !----------------------------------------------------------------------------------------------------------------------------
743    t1 => sections(is)%trac(:)
744    !--------------------------------------------------------------------------------------------------------------------------
745    DO iq=1, nt(is)                                                  !=== LOOP ON TRACERS
746    !--------------------------------------------------------------------------------------------------------------------------
747      tnam = TRIM(t1(iq)%name)                                       !--- Original name
748      IF(COUNT(t1%name == tnam) == 1) CYCLE                          !--- Current tracer is not duplicated: finished
749      tnam_new = TRIM(tnam)//'_'//TRIM(sections(is)%name)            !--- Same with section extension
750      nq = SUM(nt(1:is-1))                                           !--- Number of tracers in previous sections
751      ns = nt(is)                                                    !--- Number of tracers in the current section
752      tr(iq + nq)%name = TRIM(tnam_new)                              !--- Modify tracer name
753      WHERE(tr(1+nq:ns+nq)%parent==tnam) tr(1+nq:ns+nq)%parent=tnam_new  !--- Modify parent name
754    !--------------------------------------------------------------------------------------------------------------------------
755    END DO
756  !----------------------------------------------------------------------------------------------------------------------------
757  END DO
758  !----------------------------------------------------------------------------------------------------------------------------
759  CALL sortTracers(tr)
760END FUNCTION cumulTracers
761!==============================================================================================================================
762
763!==============================================================================================================================
764SUBROUTINE setDirectKeys(tr)
765  TYPE(trac_type), INTENT(INOUT) :: tr(:)
766
767  !--- Update %iqParent, %iqDescen, %nqDescen, %nqChilds
768  CALL indexUpdate(tr)
769
770  !--- Extract some direct-access keys
771!  DO iq = 1, SIZE(tr)
772!    tr(iq)%keys%<key> = getKey_prv(it, "<key>", tr%keys, <default_value> )
773!  END DO
774END SUBROUTINE setDirectKeys
775!==============================================================================================================================
776
777!==============================================================================================================================
778LOGICAL FUNCTION dispTraSection(message, sname, modname) RESULT(lerr)
779  CHARACTER(LEN=*), INTENT(IN) :: message, sname, modname
780  INTEGER :: idb, iq, nq
781  INTEGER, ALLOCATABLE :: hadv(:), vadv(:)
782  CHARACTER(LEN=maxlen), ALLOCATABLE :: phas(:)
783  TYPE(trac_type), POINTER :: tm(:)
784  lerr = .FALSE.
785  idb = strIdx(dBase(:)%name, sname); IF(idb == 0) RETURN
786  tm => dBase(idb)%trac
787  nq = SIZE(tm)
788  !--- BEWARE ! Can't use the "getKeyByName" functions yet.
789  !             Names must first include the phases for tracers defined on multiple lines.
790  hadv = str2int([(fgetKey(iq, 'hadv',  tm(:)%keys, '10'), iq=1, nq)])
791  vadv = str2int([(fgetKey(iq, 'vadv',  tm(:)%keys, '10'), iq=1, nq)])
792  phas =         [(fgetKey(iq, 'phases',tm(:)%keys, 'g' ), iq=1, nq)]
793  CALL msg(TRIM(message)//':', modname)
794  IF(ALL(tm(:)%parent == '')) THEN
795    IF(test(dispTable('iiiss',   ['iq    ','hadv  ','vadv  ','name  ','phase '], cat(tm%name, phas), &
796                 cat([(iq, iq=1, nq)], hadv, vadv),                 nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
797  ELSE
798    IF(test(dispTable('iiissis', ['iq    ','hadv  ','vadv  ','name  ','parent','igen  ','phase '], cat(tm%name, tm%parent, &
799      tm%phase), cat([(iq, iq=1, nq)], hadv, vadv, tm%iGeneration), nColMax=maxTableWidth, nHead=2, sub=modname), lerr)) RETURN
800  END IF
801END FUNCTION dispTraSection
802!==============================================================================================================================
803!==============================================================================================================================
804
805
806!==============================================================================================================================
807!== CREATE A SCALAR ALIAS OF THE COMPONENT OF THE TRACERS DESCRIPTOR "t" NAMED "tname" ========================================
808!==============================================================================================================================
809FUNCTION aliasTracer(tname, t) RESULT(out)
810  TYPE(trac_type),         POINTER    :: out
811  CHARACTER(LEN=*),        INTENT(IN) :: tname
812  TYPE(trac_type), TARGET, INTENT(IN) :: t(:)
813  INTEGER :: it
814  it = strIdx(t(:)%name, tname)
815  out => NULL(); IF(it /= 0) out => t(it)
816END FUNCTION aliasTracer
817!------------------------------------------------------------------------------------------------------------------------------
818
819
820!==============================================================================================================================
821!=== FROM A LIST OF INDEXES OR NAMES, CREATE A SUBSET OF THE TRACERS DESCRIPTORS LIST "trac" ==================================
822!==============================================================================================================================
823FUNCTION trSubset_Indx(trac,idx) RESULT(out)
824  TYPE(trac_type), ALLOCATABLE             ::  out(:)
825  TYPE(trac_type), ALLOCATABLE, INTENT(IN) :: trac(:)
826  INTEGER,                      INTENT(IN) ::  idx(:)
827  out = trac(idx)
828  CALL indexUpdate(out)
829END FUNCTION trSubset_Indx
830!------------------------------------------------------------------------------------------------------------------------------
831FUNCTION trSubset_Name(trac,nam) RESULT(out)
832  TYPE(trac_type), ALLOCATABLE             ::  out(:)
833  TYPE(trac_type), ALLOCATABLE, INTENT(IN) :: trac(:)
834  CHARACTER(LEN=*),             INTENT(IN) ::  nam(:)
835  out = trac(strIdx(trac(:)%name, nam))
836  CALL indexUpdate(out)
837END FUNCTION trSubset_Name
838!------------------------------------------------------------------------------------------------------------------------------
839
840
841!==============================================================================================================================
842!=== CREATE THE SUBSET OF THE TRACERS DESCRIPTORS LIST "trac" HAVING THE FIRST GENERATION ANCESTOR NAMED "nam" ================
843!==============================================================================================================================
844FUNCTION trSubset_gen0Name(trac,nam) RESULT(out)
845  TYPE(trac_type), ALLOCATABLE             ::  out(:)
846  TYPE(trac_type), ALLOCATABLE, INTENT(IN) :: trac(:)
847  CHARACTER(LEN=*),             INTENT(IN) ::  nam
848  out = trac(strFind(delPhase(trac(:)%gen0Name), nam))
849  CALL indexUpdate(out)
850END FUNCTION trSubset_gen0Name
851!------------------------------------------------------------------------------------------------------------------------------
852
853
854!==============================================================================================================================
855!=== UPDATE THE INDEXES iqParent, iqDescend AND iGeneration IN THE TRACERS DESCRIPTOR LIST "tr" (USEFULL FOR SUBSETS) =========
856!==============================================================================================================================
857SUBROUTINE indexUpdate(tr)
858  TYPE(trac_type), INTENT(INOUT) :: tr(:)
859  INTEGER :: iq, ig, ng, igen, ngen
860  INTEGER, ALLOCATABLE :: ix(:)
861  tr(:)%iqParent = strIdx( tr(:)%name, tr(:)%parent )                !--- Parent index
862  ngen = MAXVAL(tr(:)%iGeneration, MASK=.TRUE.)
863  DO iq = 1, SIZE(tr)
864    ig = tr(iq)%iGeneration
865    IF(ALLOCATED(tr(iq)%iqDescen)) DEALLOCATE(tr(iq)%iqDescen)
866    ALLOCATE(tr(iq)%iqDescen(0))
867    ix = idxAncestor(tr, igen=ig)                                    !--- Ancestor of generation "ng" for each tr
868    DO igen = ig+1, ngen
869      tr(iq)%iqDescen = [tr(iq)%iqDescen, find(ix==iq .AND. tr%iGeneration==igen)]
870      tr(iq)%nqDescen = SIZE(tr(iq)%iqDescen)
871      IF(igen == ig+1) tr(iq)%nqChilds=tr(iq)%nqDescen
872    END DO
873  END DO
874END SUBROUTINE indexUpdate
875!------------------------------------------------------------------------------------------------------------------------------
876 
877 
878!==============================================================================================================================
879!=== READ FILE "fnam" TO APPEND THE "dBase" TRACERS DATABASE WITH AS MUCH SECTIONS AS PARENTS NAMES IN "isot(:)%parent":   ====
880!===  * Each section dBase(i)%name contains the isotopes "dBase(i)%trac(:)" descending on "dBase(i)%name"="iso(i)%parent"  ====
881!===  * For each isotopes class, the <key>=<val> vector of each tracer is moved into the isotopes descriptor "isot"        ====
882!=== NOTES:                                                                                                                ====
883!===  * Most of the "isot" components have been defined in the calling routine (initIsotopes):                             ====
884!===      parent,  nzone, zone(:),  niso, keys(:)%name,  ntiso, trac(:),  nphas, phas,  iqIsoPha(:,:),  itZonPhi(:,:)      ====
885!===  * Same syntax for isotopes file and "tracer.def": a tracers section contains one line for each of its isotopes       ====
886!===  * Each tracers section can contain a "params" virtual isotope line of isotopes parameters default values             ====
887!===  * In case keys are found both in the "params" section and the "*.def" file, the later value is retained              ====
888!===  * On each isotope line, defined keys can be used for other keys defintions (single level depth substitution)         ====
889!===  * The routine gives an error if a required isotope is not available in the database stored in "fnam"                 ====
890!==============================================================================================================================
891LOGICAL FUNCTION readIsotopesFile(fnam, isot) RESULT(lerr)
892  CHARACTER(LEN=*),        INTENT(IN)    :: fnam                     !--- Input file name
893  TYPE(isot_type), TARGET, INTENT(INOUT) :: isot(:)                  !--- Isotopes descriptors (field "prnt" must be defined !)
894  INTEGER :: ik, is, it, idb, nk0, i, iis
895  INTEGER :: nk, ns, nt, ndb, nb0, i0
896  CHARACTER(LEN=maxlen), POINTER     :: k(:), v(:), k0(:), v0(:)
897  CHARACTER(LEN=maxlen), ALLOCATABLE :: vals(:)
898  CHARACTER(LEN=maxlen)              :: val, modname
899  TYPE(keys_type),           POINTER ::   ky(:)
900  TYPE(trac_type),           POINTER ::   tt(:), t
901  TYPE(dataBase_type),   ALLOCATABLE ::  tdb(:)
902  LOGICAL,               ALLOCATABLE :: liso(:)
903  modname = 'readIsotopesFile'
904
905  !--- THE INPUT FILE MUST BE PRESENT
906  IF(test(fmsg('Missing isotopes parameters file "'//TRIM(fnam)//'"', modname, testFile(fnam)),lerr)) RETURN
907
908  !--- READ THE FILE SECTIONS, ONE EACH PARENT TRACER
909  nb0 = SIZE(dBase, DIM=1)+1                                         !--- Next database element index
910  IF(test(readSections(fnam,strStack(isot(:)%parent,'|')),lerr)) RETURN !--- Read sections, one each parent tracer
911  ndb = SIZE(dBase, DIM=1)                                           !--- Current database size
912  DO idb = nb0, ndb
913   iis = idb-nb0+1
914
915    !--- GET FEW GLOBAL KEYS FROM "def" FILES AND ADD THEM TO THE 'params' SECTION
916    CALL addKeysFromDef(dBase(idb)%trac, 'params')
917
918    !--- SUBSTITUTE THE KEYS DEFINED IN THE 'params' VIRTUAL TRACER ; SUBSTITUTE LOCAL KEYS ; REMOVE 'params' VIRTUAL TRACER
919    CALL subDefault(dBase(idb)%trac, 'params', .TRUE.)
920
921    tt => dBase(idb)%trac
922
923    !--- REDUCE THE EXPRESSIONS TO OBTAIN SCALARS AND TRANSFER THEM TO THE "isot" ISOTOPES DESCRIPTORS VECTOR
924    DO it = 1, SIZE(dBase(idb)%trac)
925      t => dBase(idb)%trac(it)
926      is = strIdx(isot(iis)%keys(:)%name, t%name)                    !--- Index in "isot(iis)%keys(:)%name" of isotope "t%name"
927      IF(is == 0) CYCLE
928      liso = reduceExpr(t%keys%val, vals)                            !--- Reduce expressions (for substituted variables)
929      IF(test(ANY(liso), lerr)) RETURN                               !--- Some non-numerical elements were found
930      isot(iis)%keys(is)%key = PACK(t%keys%key, MASK=.NOT.liso)
931      isot(iis)%keys(is)%val = PACK(  vals,     MASK=.NOT.liso)
932    END DO
933
934    !--- CHECK FOR MISSING ISOTOPES (NO KEYS ALLOCATED)
935    liso = [( ALLOCATED(isot(iis)%keys(is)%key), is=1, SIZE(isot(iis)%keys) )]
936    IF(test(checkList(isot(iis)%keys(:)%name, .NOT.liso, &
937      'Check file "'//TRIM(fnam)//'" in section "'//TRIM(dBase(idb)%name)//'"', 'isotopes', 'missing'),lerr)) RETURN
938  END DO
939
940  !--- CLEAN THE DATABASE ENTRIES
941  IF(nb0 == 1) THEN
942    DEALLOCATE(dBase); ALLOCATE(dBase(0))
943  ELSE
944    ALLOCATE(tdb(nb0-1)); tdb(1:nb0-1)=dBase(1:nb0-1); CALL MOVE_ALLOC(FROM=tdb, TO=dBase)
945  END IF
946
947  !--- GET THE isoCheck ENTRY FROM THE *.DEF FILES (MIGHT BE CHANGED TO A CLASS-DEPENDANT KEYWORD)
948  CALL get_in('ok_iso_verif', isot(strIdx(isot%parent, 'H2O'))%check, .FALSE.)
949
950  lerr = dispIsotopes(isot, 'Isotopes parameters read from file "'//TRIM(fnam)//'"', modname)
951
952END FUNCTION readIsotopesFile
953!==============================================================================================================================
954
955!==============================================================================================================================
956!=== IF ISOTOPES (2ND GENERATION TRACERS) ARE DETECTED:                                                                     ===
957!===    * COMPUTE MOST OF THE RELATED QUANTITIES ("isot" COMPONENTS).                                                       ===
958!===    * COMPUTE FEW ISOTOPES-DEDICATED "trac" COMPONENTS                                                                  ===
959!===    * CALL readIsotopesFile TO GET PHYSICAL QUANTITIES (<key>=<val> PAIRS)                                              ===
960!===      NOTE: THIS IS DONE HERE (IN A ROUTINE CALLED BY THE DYNAMIC), BECAUSE THE DYNAMIC NEEDS FEW PHYSICAL PARAMETERS.  ===
961!==============================================================================================================================
962LOGICAL FUNCTION initIsotopes(trac, isot) RESULT(lerr)
963  TYPE(trac_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: trac(:)
964  TYPE(isot_type), ALLOCATABLE, TARGET, INTENT(INOUT) :: isot(:)
965  CHARACTER(LEN=maxlen), ALLOCATABLE :: p(:), str(:)                 !--- Temporary storage
966  CHARACTER(LEN=maxlen) :: iname
967  CHARACTER(LEN=1)   :: ph                                           !--- Phase
968  INTEGER :: nbIso, ic, ip, iq, it, iz
969  LOGICAL, ALLOCATABLE :: ll(:)                                      !--- Mask
970  TYPE(trac_type), POINTER   ::  t(:), t1
971  TYPE(isot_type), POINTER   ::  i
972  lerr = .FALSE.
973
974  t => trac
975
976  p = PACK(delPhase(t%parent), MASK = t%type=='tracer' .AND. t%iGeneration==1) !--- Parents of generation 1 isotopes
977  CALL strReduce(p, nbIso)
978  ALLOCATE(isot(nbIso))
979
980  IF(nbIso==0) RETURN                                                !=== NO ISOTOPES: FINISHED
981
982  !--- ISOTOPES RELATED VARIABLES ; NULL OR EMPTY IF NO ISOTOPES
983  isot(:)%parent = p
984  DO ic = 1, SIZE(p)                                                 !--- Loop on isotopes classes
985    i => isot(ic)
986    iname = i%parent                                                 !--- Current isotopes class name (parent tracer name)
987
988    !=== Isotopes childs of tracer "iname": mask, names, number (same for each phase of "iname")
989    ll = t(:)%type=='tracer' .AND. delPhase(t(:)%parent) == iname .AND. t(:)%phase == 'g'
990    str = PACK(delPhase(t(:)%name), MASK = ll)                       !--- Effectively found isotopes of "iname"
991    i%niso = SIZE(str)                                               !--- Number of "effectively found isotopes of "iname"
992    ALLOCATE(i%keys(i%niso))
993    FORALL(it = 1:i%niso) i%keys(it)%name = str(it)
994
995    !=== Geographic tagging tracers descending on tracer "iname": mask, names, number
996    ll = t(:)%type=='tag'    .AND. delPhase(t(:)%gen0Name) == iname .AND. t(:)%iGeneration == 2
997    i%zone = PACK(strTail(t(:)%name,'_'), MASK = ll)                 !--- Tagging zones names  for isotopes category "iname"
998    CALL strReduce(i%zone)
999    i%nzone = SIZE(i%zone)                                           !--- Tagging zones number for isotopes category "iname"
1000
1001    !=== Geographic tracers of the isotopes childs of tracer "iname" (same for each phase of "iname")
1002    !    NOTE: One might like to create a similar variable for 2nd generation tagging tracers (tagging the gen1 tracers)
1003    str = PACK(delPhase(t(:)%name), MASK=ll)
1004    CALL strReduce(str)
1005    i%ntiso = i%niso + SIZE(str)                                     !--- Number of isotopes + their geographic tracers [ntiso]
1006    ALLOCATE(i%trac(i%ntiso))
1007    FORALL(it = 1:i%niso) i%trac(it) = i%keys(it)%name
1008    FORALL(it = i%niso+1:i%ntiso) i%trac(it) = str(it-i%niso)
1009
1010    !=== Phases for tracer "iname"
1011    i%phase = ''
1012    DO ip = 1, nphases; ph = known_phases(ip:ip); IF(strIdx(t%name,addPhase(iname, ph)) /= 0) i%phase = TRIM(i%phase)//ph; END DO
1013    i%nphas = LEN_TRIM(i%phase)                                       !--- Equal to "nqo" for water
1014
1015    !=== Tables giving the index in a table of effectively found items for each dynamical tracer (1<=iq<=nqtot)
1016    DO iq = 1, SIZE(t)
1017      t1 => trac(iq)
1018      IF(delPhase(t1%gen0Name)/=iname .OR. t1%iGeneration==0) CYCLE  !--- Only deal with tracers descending on "iname"
1019      t1%iso_iGroup = ic                                             !--- Isotopes family       idx in list "isotopes(:)%parent"
1020      t1%iso_iName  = strIdx(i%trac, delPhase(strHead(t1%name,'_'))) !--- Current isotope       idx in effective isotopes list
1021      t1%iso_iZone  = strIdx(i%zone,          strTail(t1%name,'_') ) !--- Current isotope zone  idx in effective zones    list
1022      t1%iso_iPhase =  INDEX(i%phase,TRIM(t1%phase))                 !--- Current isotope phase idx in effective phases   list
1023      IF(t1%iGeneration /= 2) t1%iso_iZone = 0                       !--- Skip possible generation 1 tagging tracers
1024    END DO
1025
1026    !=== Table used to get iq (index in dyn array, size nqtot) from the isotope and phase indexes ; the full isotopes list
1027    !    (including tagging tracers) is sorted this way:  iso1, iso2, ..., iso1_zone1, iso2_zone1, ..., iso1_zoneN, iso2_zoneN
1028    i%iqIsoPha = RESHAPE( [( (strIdx(t%name,  addPhase(i%trac(it),i%phase(ip:ip))),    it=1, i%ntiso), ip=1, i%nphas)], &
1029                         [i%ntiso, i%nphas] )
1030    !=== Table used to get ix (index in tagging tracers isotopes list, size ntiso) from the zone and isotope indexes
1031    i%itZonIso = RESHAPE( [( (strIdx(i%trac(:), TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))), iz=1, i%nzone), it=1, i%niso )], &
1032                         [i%nzone, i%niso] )
1033  END DO
1034
1035  !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE
1036  !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat AND alpha_ideal)
1037  lerr = readIsotopesFile('isotopes_params.def',isot)
1038
1039END FUNCTION initIsotopes
1040!==============================================================================================================================
1041
1042
1043!==============================================================================================================================
1044LOGICAL FUNCTION dispIsotopes(ides, message, modname) RESULT(lerr)
1045  TYPE(isot_type),  INTENT(IN) :: ides(:)                            !--- Isotopes descriptor vector
1046  CHARACTER(LEN=*), INTENT(IN) :: message                            !--- Message to display
1047  CHARACTER(LEN=*), INTENT(IN) :: modname                            !--- Calling subroutine name
1048  INTEGER :: ik, nk, ip, it, nt
1049  CHARACTER(LEN=maxlen) :: prf
1050  CHARACTER(LEN=maxlen), ALLOCATABLE :: ttl(:), val(:,:)
1051  CALL msg(TRIM(message)//':', modname)
1052  DO ip = 1, SIZE(ides)                                              !--- Loop on parents tracers
1053    nk = SIZE(ides(ip)%keys(1)%key)                                  !--- Same keys for each isotope
1054    nt = SIZE(ides(ip)%keys)                                         !--- Number of isotopes
1055    prf = 'i'//REPEAT('s',nk+1)                                      !--- Profile for table printing
1056    ALLOCATE(ttl(nk+2), val(nt,nk+1))
1057    ttl(1:2) = ['iq  ','name']; ttl(3:nk+2) = ides(ip)%keys(1)%key(:)!--- Titles line with keys names
1058    val(:,1) = ides(ip)%keys(:)%name                                 !--- Values table 1st column: isotopes names 
1059    DO ik = 1, nk
1060      DO it = 1, nt
1061        val(it,ik+1) = ides(ip)%keys(it)%val(ik)                     !--- Other columns: keys values
1062      END DO
1063    END DO
1064    IF(test(fmsg('Problem with the table content', modname, dispTable(prf, ttl, val, &
1065            cat([(it,it=1,nt)]), rFmt='(EN8.4)', nColMax=maxTableWidth, nHead=2, sub=modname)), lerr)) RETURN
1066    DEALLOCATE(ttl, val)
1067  END DO       
1068END FUNCTION dispIsotopes
1069!==============================================================================================================================
1070
1071
1072!==============================================================================================================================
1073SUBROUTINE addKey_1(key, val, ky, lOverWrite)
1074!------------------------------------------------------------------------------------------------------------------------------
1075! Purpose: Add the <key>=<val> pair in the "ky" keys descriptor.
1076!------------------------------------------------------------------------------------------------------------------------------
1077  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
1078  TYPE(keys_type),   INTENT(INOUT) :: ky
1079  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
1080  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
1081  INTEGER :: iky, nky
1082  LOGICAL :: lo
1083!------------------------------------------------------------------------------------------------------------------------------
1084  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
1085  iky = strIdx(ky%key,key)
1086  IF(iky == 0) THEN
1087    nky = SIZE(ky%key)
1088    IF(nky == 0) THEN; ky%key = [key]; ky%val = [val]; ELSE; ky%key = [ky%key, key]; ky%val = [ky%val, val]; END IF
1089  ELSE IF(lo) THEN                                                   !--- Overwriting
1090    ky%key(iky) = key; ky%val(iky) = val
1091  END IF
1092END SUBROUTINE addKey_1
1093!==============================================================================================================================
1094SUBROUTINE addKey_m(key, val, ky, lOverWrite)
1095!------------------------------------------------------------------------------------------------------------------------------
1096! Purpose: Add the <key>=<val> pair in all the components of the "ky" keys descriptor.
1097!------------------------------------------------------------------------------------------------------------------------------
1098  CHARACTER(LEN=*),  INTENT(IN)    :: key, val
1099  TYPE(keys_type),   INTENT(INOUT) :: ky(:)
1100  LOGICAL, OPTIONAL, INTENT(IN)    :: lOverWrite
1101  INTEGER :: itr
1102  LOGICAL :: lo
1103!------------------------------------------------------------------------------------------------------------------------------
1104  lo=.FALSE.; IF(PRESENT(lOverWrite)) lo=lOverWrite
1105  DO itr = 1, SIZE(ky); CALL addKey_1(key, val, ky(itr), lo); END DO
1106END SUBROUTINE addKey_m
1107!==============================================================================================================================
1108SUBROUTINE addKeysFromDef(t, tr0)
1109!------------------------------------------------------------------------------------------------------------------------------
1110! Purpose: The values of the keys of the tracer named "tr0" are overwritten by the values found in the *.def files, if any.
1111!------------------------------------------------------------------------------------------------------------------------------
1112  TYPE(trac_type), ALLOCATABLE, INTENT(INOUT) :: t(:)
1113  CHARACTER(LEN=*),             INTENT(IN)    :: tr0
1114  CHARACTER(LEN=maxlen) :: val
1115  INTEGER               :: ik, jd
1116  jd = strIdx(t%name, tr0)
1117  IF(jd == 0) RETURN
1118  DO ik = 1, SIZE(t(jd)%keys%key)
1119    CALL get_in(t(jd)%keys%key(ik), val, '*none*')
1120    IF(val /= '*none*') CALL addKey_1(t(jd)%keys%key(ik), val, t(jd)%keys, .TRUE.)
1121  END DO
1122END SUBROUTINE addKeysFromDef
1123!==============================================================================================================================
1124SUBROUTINE delKey_1(itr, keyn, ky)
1125!------------------------------------------------------------------------------------------------------------------------------
1126! Purpose: Internal routine.
1127!   Remove <key>=<val> pairs in the "itr"th component of the "ky" keys descriptor.
1128!------------------------------------------------------------------------------------------------------------------------------
1129  INTEGER,          INTENT(IN)    :: itr
1130  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
1131  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
1132  CHARACTER(LEN=maxlen), ALLOCATABLE :: k(:), v(:)
1133  LOGICAL,               ALLOCATABLE :: ll(:)
1134  INTEGER :: iky
1135!------------------------------------------------------------------------------------------------------------------------------
1136  IF(itr<=0 .OR. itr>SIZE(ky, DIM=1)) RETURN                          !--- Index is out of range
1137  ll = [( ALL(keyn/=ky(itr)%keys%key(iky)), iky=1, SIZE(ky(itr)%keys%key) )]
1138  k = PACK(ky(itr)%keys%key, MASK=ll); CALL MOVE_ALLOC(FROM=k, TO=ky(itr)%keys%key)
1139  v = PACK(ky(itr)%keys%val, MASK=ll); CALL MOVE_ALLOC(FROM=v, TO=ky(itr)%keys%val)
1140END SUBROUTINE delKey_1
1141!==============================================================================================================================
1142SUBROUTINE delKey(keyn, ky)
1143!------------------------------------------------------------------------------------------------------------------------------
1144! Purpose: Internal routine.
1145!   Remove <key>=<val> pairs in all the components of the "t" tracers descriptor.
1146!------------------------------------------------------------------------------------------------------------------------------
1147  CHARACTER(LEN=*), INTENT(IN)    :: keyn(:)
1148  TYPE(trac_type),  INTENT(INOUT) :: ky(:)
1149  INTEGER :: iky
1150!------------------------------------------------------------------------------------------------------------------------------
1151  DO iky = 1, SIZE(ky); CALL delKey_1(iky, keyn, ky); END DO
1152END SUBROUTINE delKey
1153!==============================================================================================================================
1154
1155
1156!==============================================================================================================================
1157!=== PUBLIC ROUTINES: GET A KEY FROM A <key>=<val> LIST ; VECTORS, TRACER AND DATABASE VERSIONS ===============================
1158!=== BEWARE !!! IF THE "ky" ARGUMENT IS NOT PRESENT, THEN THE VARIABLES "tracers" AND "isotopes" ARE USED. ====================
1159!===     THEY ARE LOCAL TO THIS MODULE, SO MUST MUST BE INITIALIZED FIRST USING the "getKey_init" ROUTINE  ====================
1160!==============================================================================================================================
1161SUBROUTINE getKey_init(tracers_, isotopes_)
1162  TYPE(trac_type), OPTIONAL, INTENT(IN) ::  tracers_(:)
1163  TYPE(isot_type), OPTIONAL, INTENT(IN) :: isotopes_(:)
1164  IF(PRESENT( tracers_))  tracers =  tracers_
1165  IF(PRESENT(isotopes_)) isotopes = isotopes_
1166END SUBROUTINE getKey_init
1167!==============================================================================================================================
1168CHARACTER(LEN=maxlen) FUNCTION fgetKeyByIndex_s1(itr, keyn, ky, def_val) RESULT(val)
1169!------------------------------------------------------------------------------------------------------------------------------
1170! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer index.
1171!------------------------------------------------------------------------------------------------------------------------------
1172  INTEGER,                    INTENT(IN) :: itr
1173  CHARACTER(LEN=*),           INTENT(IN) :: keyn
1174  TYPE(keys_type),            INTENT(IN) :: ky(:)
1175  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: def_val
1176!------------------------------------------------------------------------------------------------------------------------------
1177  INTEGER :: iky
1178  iky = 0;  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
1179  val = ''; IF(iky /= 0) val = ky(itr)%val(iky)                      !--- Key was found
1180  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
1181END FUNCTION fgetKeyByIndex_s1
1182!==============================================================================================================================
1183CHARACTER(LEN=maxlen) FUNCTION fgetKeyByName_s1(tname, keyn, ky, def_val, lerr) RESULT(val)
1184!------------------------------------------------------------------------------------------------------------------------------
1185! Purpose: Internal function ; get a string formatted key value (returned argument) from its name and the tracer name.
1186!------------------------------------------------------------------------------------------------------------------------------
1187  CHARACTER(LEN=*),           INTENT(IN)  :: tname, keyn
1188  TYPE(keys_type),            INTENT(IN)  :: ky(:)
1189  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: def_val
1190  LOGICAL,          OPTIONAL, INTENT(OUT) :: lerr
1191!------------------------------------------------------------------------------------------------------------------------------
1192  INTEGER :: iky, itr
1193  val = ''; iky = 0
1194  itr = strIdx(ky(:)%name, tname)                                    !--- Get the index of the wanted tracer
1195  IF(PRESENT(lerr)) lerr = itr==0; IF(itr == 0) RETURN
1196  IF(itr >  0 .AND. itr <= SIZE(ky)) iky = strIdx(ky(itr)%key(:), keyn)
1197  IF(iky /= 0) val = ky(itr)%val(iky)                                !--- Key was found
1198  IF(PRESENT(def_val) .AND. iky == 0) val = def_val                  !--- Default value from arguments
1199END FUNCTION fgetKeyByName_s1
1200!==============================================================================================================================
1201LOGICAL FUNCTION getKeyByName_s1(keyn, val, tname, ky) RESULT(lerr)
1202  !--- Purpose: Get the value of the key named "keyn" for the tracer named "tnam".
1203  !     * "ky" unspecified: try in "tracers" for "tnam" with phase and tagging suffixes, then in "isotopes" without.
1204  !     * "ky"   specified: try in "ky"      for "tnam" with phase and tagging suffixes, then without.
1205  !    The returned error code is always .FALSE.: an empty string is returned when the key hasn't been found.
1206  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
1207  CHARACTER(LEN=maxlen),     INTENT(OUT) :: val
1208  CHARACTER(LEN=*),          INTENT(IN)  :: tname
1209  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
1210  CHARACTER(LEN=maxlen) :: tnam
1211  INTEGER, ALLOCATABLE  :: is(:)
1212  INTEGER :: i, itr
1213  tnam = delPhase(strHead(tname,'_',.FALSE.))                        !--- Remove tag and phase
1214  IF(PRESENT(ky)) THEN
1215    val = fgetKeyByName_s1(tname, keyn, ky, lerr=lerr)               !--- "ky" and "tname"
1216    IF(val /= '' .OR. lerr)      RETURN
1217    val = fgetKeyByName_s1(tnam,  keyn, ky, lerr=lerr)               !--- "ky" and "tnam"
1218  ELSE
1219    IF(.NOT.ALLOCATED(tracers))  RETURN
1220    val = fgetKeyByName_s1(tname, keyn, tracers(:)%keys, lerr=lerr)  !--- "tracers" and "tname"
1221    IF(val /= ''.AND..NOT.lerr)  RETURN
1222    IF(.NOT.ALLOCATED(isotopes)) RETURN
1223    IF(SIZE(isotopes) == 0)      RETURN
1224    !--- Search the "is" isotopes class index of the isotope named "tnam"
1225    is = find([(ANY(isotopes(i)%keys(:)%name == tnam), i=1, SIZE(isotopes))])
1226    IF(test(SIZE(is) == 0,lerr)) RETURN
1227    val = fgetKeyByName_s1(tname, keyn, isotopes(is(1))%keys(:),lerr=lerr)!--- "isotopes" and "tnam"
1228  END IF
1229END FUNCTION getKeyByName_s1
1230!==============================================================================================================================
1231LOGICAL FUNCTION getKeyByName_sm(keyn, val, tname, ky) RESULT(lerr)
1232  CHARACTER(LEN=*),                   INTENT(IN)  :: keyn
1233  CHARACTER(LEN=maxlen), ALLOCATABLE, INTENT(OUT) ::   val(:)
1234  CHARACTER(LEN=*),         OPTIONAL, INTENT(IN)  :: tname(:)
1235  TYPE(keys_type),          OPTIONAL, INTENT(IN)  ::    ky(:)
1236  TYPE(keys_type),           POINTER :: k(:)
1237  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
1238  INTEGER :: iq, nq
1239  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
1240  IF(PRESENT(ky   )) nq = SIZE(ky%name)
1241  IF(PRESENT(tname)) nq = SIZE(  tname)
1242  ALLOCATE(val(nq))
1243  IF(PRESENT(tname)) THEN
1244    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
1245    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_s1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
1246  ELSE;                  lerr = ANY([(getKeyByName_s1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
1247  END IF
1248END FUNCTION getKeyByName_sm
1249!==============================================================================================================================
1250LOGICAL FUNCTION getKeyByName_i1(keyn, val, tname, ky) RESULT(lerr)
1251  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
1252  INTEGER,                   INTENT(OUT) :: val
1253  CHARACTER(LEN=*),          INTENT(IN)  :: tname
1254  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
1255  CHARACTER(LEN=maxlen) :: sval
1256  INTEGER :: ierr
1257  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
1258  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
1259  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',        modname, lerr), lerr)) RETURN
1260  READ(sval, *, IOSTAT=ierr) val
1261  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not an integer', modname, lerr), lerr)) RETURN
1262END FUNCTION getKeyByName_i1
1263!==============================================================================================================================
1264LOGICAL FUNCTION getKeyByName_im(keyn, val, tname, ky) RESULT(lerr)
1265  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
1266  INTEGER,       ALLOCATABLE, INTENT(OUT) ::   val(:)
1267  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: tname(:)
1268  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::    ky(:)
1269  TYPE(keys_type),           POINTER :: k(:)
1270  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
1271  INTEGER :: iq, nq
1272  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
1273  IF(PRESENT(ky   )) nq = SIZE(ky%name)
1274  IF(PRESENT(tname)) nq = SIZE(  tname)
1275  ALLOCATE(val(nq))
1276  IF(PRESENT(tname)) THEN
1277    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
1278    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_i1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
1279  ELSE;                  lerr = ANY([(getKeyByName_i1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
1280  END IF
1281END FUNCTION getKeyByName_im
1282!==============================================================================================================================
1283LOGICAL FUNCTION getKeyByName_r1(keyn, val, tname, ky) RESULT(lerr)
1284  CHARACTER(LEN=*),          INTENT(IN)  :: keyn
1285  REAL,                      INTENT(OUT) :: val
1286  CHARACTER(LEN=*),          INTENT(IN)  :: tname
1287  TYPE(keys_type), OPTIONAL, INTENT(IN)  :: ky(:)
1288  CHARACTER(LEN=maxlen) :: sval
1289  INTEGER :: ierr
1290  IF(     PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname, ky)
1291  IF(.NOT.PRESENT(ky)) lerr = getKeyByName_s1(keyn, sval, tname)
1292  IF(test(fmsg('key "'//TRIM(keyn)//'" or tracer "'//TRIM(tname)//'" is missing',    modname, lerr), lerr)) RETURN
1293  READ(sval, *, IOSTAT=ierr) val
1294  IF(test(fmsg('key "'//TRIM(keyn)//'" of tracer "'//TRIM(tname)//'" is not a real', modname, lerr), lerr)) RETURN
1295END FUNCTION getKeyByName_r1
1296!==============================================================================================================================
1297LOGICAL FUNCTION getKeyByName_rm(keyn, val, tname, ky) RESULT(lerr)
1298  CHARACTER(LEN=*),           INTENT(IN)  :: keyn
1299  REAL,          ALLOCATABLE, INTENT(OUT) ::   val(:)
1300  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: tname(:)
1301  TYPE(keys_type),  OPTIONAL, INTENT(IN)  ::    ky(:)
1302  TYPE(keys_type),           POINTER :: k(:)
1303  CHARACTER(LEN=maxlen), ALLOCATABLE :: n(:)
1304  INTEGER :: iq, nq
1305  IF(test(.NOT.(PRESENT(tname).OR.PRESENT(ky)), lerr)) RETURN
1306  IF(PRESENT(ky   )) nq = SIZE(ky%name)
1307  IF(PRESENT(tname)) nq = SIZE(  tname)
1308  ALLOCATE(val(nq))
1309  IF(PRESENT(tname)) THEN
1310    IF(     PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq),   ky), iq=1, nq)])
1311    IF(.NOT.PRESENT(ky)) lerr = ANY([(getKeyByName_r1(keyn, val(iq), tname(iq)      ), iq=1, nq)])
1312  ELSE;                  lerr = ANY([(getKeyByName_r1(keyn, val(iq), ky(iq)%name, ky), iq=1, nq)])
1313  END IF
1314END FUNCTION getKeyByName_rm
1315!==============================================================================================================================
1316
1317
1318!==============================================================================================================================
1319!=== REMOVE, IF ANY, OR ADD THE PHASE SUFFIX OF THE TRACER NAMED "s" ==========================================================
1320!==============================================================================================================================
1321ELEMENTAL CHARACTER(LEN=maxlen) FUNCTION delPhase(s) RESULT(out)
1322  CHARACTER(LEN=*), INTENT(IN) :: s
1323  INTEGER :: l, i, ix
1324  CHARACTER(LEN=maxlen) :: sh, st
1325  out = s
1326  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
1327
1328  !--- Special case: old phases for water, no phases separator
1329  i = INDEX(s,'_'); sh = s; IF(i/=0) sh=s(1:i-1); st='H2O'; IF(i/=0) st='H2O_'//s(i+1:LEN_TRIM(s))
1330  IF(ANY([('H2O'//old_phases(ix:ix), ix=1, nphases)] == sh)) THEN; out=st; RETURN; END IF
1331
1332  !--- Index of found phase in "known_phases"
1333  ix = MAXLOC( [(i, i=1,nphases)], MASK=[( INDEX(s, phases_sep//known_phases(i:i))/=0, i=1, nphases)], DIM=1 )
1334  IF(ix == 0) RETURN                                                           !--- No phase pattern found
1335  i = INDEX(s, phases_sep//known_phases(ix:ix))                                !--- Index of <sep><pha> pattern in "str"
1336  l = LEN_TRIM(s)
1337  IF(i == l-1) THEN                                                            !--- <var><sep><pha>       => return <var>
1338    out = s(1:l-2)
1339  ELSE IF(s(i+2:i+2) == '_') THEN                                              !--- <var><sep><pha>_<tag> => return <var>_<tag>
1340    out = s(1:i-1)//s(i+2:l)
1341  END IF
1342END FUNCTION delPhase
1343!------------------------------------------------------------------------------------------------------------------------------
1344CHARACTER(LEN=maxlen) FUNCTION addPhase_s1(s,pha) RESULT(out)
1345  CHARACTER(LEN=*),           INTENT(IN) :: s
1346  CHARACTER(LEN=1),           INTENT(IN) :: pha
1347  INTEGER :: l, i
1348  out = s
1349  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
1350  i = INDEX(s, '_')                                                            !--- /=0 for <var>_<tag> tracers names
1351  l = LEN_TRIM(s)
1352  IF(i == 0) out =  TRIM(s)//phases_sep//pha                                   !--- <var>       => return <var><sep><pha>
1353  IF(i /= 0) out = s(1:i-1)//phases_sep//pha//'_'//s(i+1:l)                    !--- <var>_<tag> => return <var><sep><pha>_<tag>
1354END FUNCTION addPhase_s1
1355!------------------------------------------------------------------------------------------------------------------------------
1356FUNCTION addPhase_sm(s,pha) RESULT(out)
1357  CHARACTER(LEN=*),           INTENT(IN) :: s(:)
1358  CHARACTER(LEN=1),           INTENT(IN) :: pha
1359  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
1360  INTEGER :: k
1361  out = [( addPhase_s1(s(k), pha), k=1, SIZE(s) )]
1362END FUNCTION addPhase_sm
1363!------------------------------------------------------------------------------------------------------------------------------
1364CHARACTER(LEN=maxlen) FUNCTION addPhase_i1(s,ipha,phases) RESULT(out)
1365  CHARACTER(LEN=*),           INTENT(IN) :: s
1366  INTEGER,                    INTENT(IN) :: ipha
1367  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
1368  out = s
1369  IF(s == '') RETURN                                                           !--- Empty string: nothing to do
1370  IF(ipha==0) RETURN                                                           !--- Null index: no phase to add
1371  IF(     PRESENT(phases)) out = addPhase_s1(s,       phases(ipha:ipha))
1372  IF(.NOT.PRESENT(phases)) out = addPhase_s1(s, known_phases(ipha:ipha))
1373END FUNCTION addPhase_i1
1374!------------------------------------------------------------------------------------------------------------------------------
1375FUNCTION addPhase_im(s,ipha,phases) RESULT(out)
1376  CHARACTER(LEN=*),           INTENT(IN) :: s(:)
1377  INTEGER,                    INTENT(IN) :: ipha
1378  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phases
1379  CHARACTER(LEN=maxlen),     ALLOCATABLE :: out(:)
1380  INTEGER :: k
1381  IF(     PRESENT(phases)) out = [( addPhase_i1(s(k), ipha,       phases), k=1, SIZE(s) )]
1382  IF(.NOT.PRESENT(phases)) out = [( addPhase_i1(s(k), ipha, known_phases), k=1, SIZE(s) )]
1383END FUNCTION addPhase_im
1384!------------------------------------------------------------------------------------------------------------------------------
1385
1386
1387!==============================================================================================================================
1388!=== GET PHASE INDEX IN THE POSSIBLE PHASES LIST OR IN A SPECIFIED LIST ("phases") ============================================
1389!==============================================================================================================================
1390INTEGER FUNCTION getiPhase(tname, phases) RESULT(iPhase)
1391  CHARACTER(LEN=*),           INTENT(IN)  :: tname
1392  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
1393  CHARACTER(LEN=maxlen) :: phase
1394  IF(     PRESENT(phases)) phase = getPhase(tname,       phases, iPhase)
1395  IF(.NOT.PRESENT(phases)) phase = getPhase(tname, known_phases, iPhase)
1396END FUNCTION getiPhase
1397!------------------------------------------------------------------------------------------------------------------------------
1398CHARACTER(LEN=1) FUNCTION getPhase(tname, phases, iPhase) RESULT(phase)
1399  CHARACTER(LEN=*),           INTENT(IN)  :: tname
1400  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: phases
1401  INTEGER,          OPTIONAL, INTENT(OUT) :: iPhase
1402  INTEGER :: ip
1403  phase = TRIM(strHead(strTail(tname, phases_sep, .TRUE.), '_', .TRUE.))
1404  IF(     PRESENT(phases)) ip = INDEX(      phases, phase)
1405  IF(.NOT.PRESENT(phases)) ip = INDEX(known_phases, phase)
1406  IF(ip == 0) phase = 'g'
1407  IF(PRESENT(iPhase)) iPhase = ip
1408END FUNCTION getPhase
1409!------------------------------------------------------------------------------------------------------------------------------
1410
1411
1412!------------------------------------------------------------------------------------------------------------------------------
1413CHARACTER(LEN=maxlen) FUNCTION old2newName_1(oldName, iPhase) RESULT(newName)
1414  !--- Convert an old style name into a new one.
1415  !    Only usable with old style "traceur.def" files, in which only water isotopes are allowed.
1416  !    In these files, H2O descendants names are: H2O<phase>[_<isotope>][_<tag>], with:
1417  !    phase = v, l or i ; isotope = eau, HDO, O18, O17 or HTO.
1418  CHARACTER(LEN=*),  INTENT(IN)  :: oldName
1419  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
1420  CHARACTER(LEN=maxlen), ALLOCATABLE :: tmp(:)
1421  INTEGER :: ix, ip, it, nt
1422  LOGICAL :: lerr, lH2O
1423  newName = oldName
1424  IF(PRESENT(iPhase)) iPhase = 1                                               !--- Default: gaseous phase
1425  lH2O=.FALSE.
1426  IF(LEN_TRIM(oldName) > 3) THEN
1427    lH2O = oldName(1:3)=='H2O' .AND. INDEX(old_phases,oldName(4:4))/=0         !--- H2O<phase>*,  with phase=="v", "l", "i" or "r"
1428    IF(LEN_TRIM(oldName) > 4) lH2O = lH2O .AND. oldName(5:5) == '_'            !--- H2O<phase>_*, with phase=="v", "l", "i" or "r"
1429  END IF
1430  IF(.NOT.lH2O) RETURN
1431  IF(LEN_TRIM(oldName)>3) THEN; IF(INDEX(old_Phases,oldName(4:4))==0) RETURN; END IF
1432  lerr = strParse(oldName, '_', tmp, n=nt)
1433  ip = strIdx([('H2O'//old_phases(ip:ip), ip=1, nphases)], tmp(1))             !--- Phase index (/=0 if any)
1434  IF(PRESENT(iPhase)) iPhase = ip
1435  newName = addPhase('H2O', ip)                                                !--- Water
1436  IF(nt == 1) RETURN                                                           !--- Water: finished
1437  ix = strIdx(oldH2OIso, tmp(2))                                               !--- Index in the known isotopes list
1438  IF(ix == 0) newName = addPhase(tmp(2),        ip)                            !--- Not an isotope
1439  IF(ix /= 0) newName = addPhase(newH2OIso(ix), ip)                            !--- Isotope
1440  IF(nt == 3) newName = TRIM(newName)//'_'//TRIM(tmp(3))                       !--- Tagging tracer
1441END FUNCTION old2newName_1
1442!------------------------------------------------------------------------------------------------------------------------------
1443FUNCTION old2newName_m(oldName, iPhase) RESULT(newName)
1444  CHARACTER(LEN=*),  INTENT(IN)  :: oldName(:)
1445  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
1446  CHARACTER(LEN=maxlen)          :: newName(SIZE(oldName))
1447  INTEGER :: i
1448  newName = [(old2newName_1(oldName(i), iPhase), i=1, SIZE(oldName))]
1449END FUNCTION old2newName_m
1450!------------------------------------------------------------------------------------------------------------------------------
1451
1452!------------------------------------------------------------------------------------------------------------------------------
1453CHARACTER(LEN=maxlen) FUNCTION new2oldName_1(newName, iPhase) RESULT(oldName)
1454  !--- Convert a new style name into an old one.
1455  !    Only convertable names are water descendants names H2O_<phase>, <isotope>_<phase>, <isotope>_<phase>_<tag>, with:
1456  !    phase = g, l or s ; isotope = H2[16]O, H[2]O, H2<[18]O, H2[17]O or H[3]O.
1457  CHARACTER(LEN=*),  INTENT(IN)    :: newName
1458  INTEGER, OPTIONAL, INTENT(OUT)   :: iPhase
1459  INTEGER :: ix, ip, it, nt
1460  LOGICAL :: lH2O
1461  CHARACTER(LEN=maxlen) :: tag
1462  ix = strIdx([(addPhase('H2O',ip), ip=1, nphases)], newName)                  !--- Phase index for H2O_<phase>
1463  IF(ix /= 0) THEN; oldName = 'H2O'//old_phases(ix:ix); RETURN; END IF         !--- H2O_<phase> case
1464  ix = strIdx(newH2OIso, strHead(newName, phases_sep, .TRUE.))                 !--- Isotope index
1465  IF(ix == 0) THEN; oldName = newName;                  RETURN; END IF         !--- Not a water descendant
1466  ip = getiPhase(newName)                                                      !--- Phase index
1467  oldName = TRIM(oldH2OIso(ix))//old_phases(ip:ip)                             !--- <isotope>_<phase>
1468  tag = strTail(delPhase(newName), TRIM(newH2OIso(ix)))                        !--- Get "_<tag>" if any
1469  IF(tag /= delPhase(newName) .AND. tag /= '') oldName = TRIM(oldName)//tag    !--- Tagging tracer
1470END FUNCTION new2oldName_1
1471!------------------------------------------------------------------------------------------------------------------------------
1472FUNCTION new2oldName_m(newName, iPhase) RESULT(oldName)
1473  CHARACTER(LEN=*),  INTENT(IN)  :: newName(:)
1474  INTEGER, OPTIONAL, INTENT(OUT) :: iPhase
1475  CHARACTER(LEN=maxlen)          :: oldName(SIZE(newName))
1476  INTEGER :: i
1477  oldName = [(new2oldName_1(newName(i), iPhase), i=1, SIZE(newName))]
1478END FUNCTION new2oldName_m
1479!------------------------------------------------------------------------------------------------------------------------------
1480
1481
1482!==============================================================================================================================
1483!=== GET THE NAME(S) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =======
1484!==============================================================================================================================
1485CHARACTER(LEN=maxlen) FUNCTION ancestor_1(t, tname, igen) RESULT(out)
1486  TYPE(trac_type),   INTENT(IN) :: t(:)
1487  CHARACTER(LEN=*),  INTENT(IN) :: tname
1488  INTEGER, OPTIONAL, INTENT(IN) :: igen
1489  INTEGER :: ig, ix
1490  ig = 0; IF(PRESENT(igen)) ig = igen
1491  ix = idxAncestor_1(t, tname, ig)
1492  out = ''; IF(ix /= 0) out = t(ix)%name
1493END FUNCTION ancestor_1
1494!------------------------------------------------------------------------------------------------------------------------------
1495FUNCTION ancestor_m(t, tname, igen) RESULT(out)
1496  CHARACTER(LEN=maxlen), ALLOCATABLE     ::   out(:)
1497  TYPE(trac_type),            INTENT(IN) ::     t(:)
1498  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
1499  INTEGER,          OPTIONAL, INTENT(IN) :: igen
1500  INTEGER, ALLOCATABLE :: ix(:)
1501  INTEGER :: ig
1502  ig = 0; IF(PRESENT(igen)) ig = igen
1503  IF(     PRESENT(tname)) ix = idxAncestor_m(t, tname,     ig)
1504  IF(.NOT.PRESENT(tname)) ix = idxAncestor_m(t, t(:)%name, ig)
1505  ALLOCATE(out(SIZE(ix))); out(:) = ''
1506  WHERE(ix /= 0) out = t(ix)%name
1507END FUNCTION ancestor_m
1508!==============================================================================================================================
1509
1510
1511!==============================================================================================================================
1512!=== GET THE INDEX(ES) OF THE ANCESTOR(S) OF TRACER(S) "tname" AT GENERATION "igen"  IN THE TRACERS DESCRIPTORS LIST "tr" =====
1513!==============================================================================================================================
1514INTEGER FUNCTION idxAncestor_1(t, tname, igen) RESULT(out)
1515! Return the name of the generation "igen" (>=0) ancestor of "tname"
1516  TYPE(trac_type),   INTENT(IN) :: t(:)
1517  CHARACTER(LEN=*),  INTENT(IN) :: tname
1518  INTEGER, OPTIONAL, INTENT(IN) :: igen
1519  INTEGER :: ig
1520  ig = 0; IF(PRESENT(igen)) ig = igen
1521  out = strIdx(t(:)%name, tname)
1522  IF(out == 0)                 RETURN            !--- Tracer not found
1523  IF(t(out)%iGeneration <= ig) RETURN            !--- Tracer has a lower generation number than asked generation 'igen"
1524  DO WHILE(t(out)%iGeneration > ig); out = strIdx(t(:)%name, t(out)%parent); END DO
1525END FUNCTION idxAncestor_1
1526!------------------------------------------------------------------------------------------------------------------------------
1527FUNCTION idxAncestor_m(t, tname, igen) RESULT(out)
1528  INTEGER,          ALLOCATABLE          ::   out(:)
1529  TYPE(trac_type),            INTENT(IN) ::     t(:)
1530  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: tname(:)
1531  INTEGER,          OPTIONAL, INTENT(IN) :: igen
1532  INTEGER :: ig, ix
1533  ig = 0; IF(PRESENT(igen)) ig = igen
1534  IF(     PRESENT(tname)) out = [(idxAncestor_1(t, tname(ix),  ig), ix=1, SIZE(tname))]
1535  IF(.NOT.PRESENT(tname)) out = [(idxAncestor_1(t, t(ix)%name, ig), ix=1, SIZE(t))]
1536END FUNCTION idxAncestor_m
1537!==============================================================================================================================
1538
1539
1540END MODULE readTracFiles_mod
Note: See TracBrowser for help on using the repository browser.