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

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

Fix in sortTracer for gfortran compatibility.

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