source: LMDZ6/trunk/libf/misc/readTracFiles_mod.f90 @ 4202

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

Remove useless "TARGET" attributes in getKeyByName_?m functions, because pointers are no longer used.

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