source: readTracFiles_mod.f90 @ 19

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

Sections names separator is no longer a coma (","), but a pipe ("|") because ioipsl "getin" is removing the comas.

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