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

Last change on this file since 4138 was 4138, checked in by acozic, 2 years ago

fixed two bug for coupling with inca
1- allow the possibility of having tracers names with H2O in it (like H2O2)
2- allow the possibility of having other water tracers in inca

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