source: readTracFiles_mod.f90 @ 6

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