source: LMDZ6/trunk/libf/phylmdiso/dyn1d/infotrac.F90 @ 3927

Last change on this file since 3927 was 3927, checked in by Laurent Fairhead, 3 years ago

Initial import of the physics wih isotopes from Camille Risi
CR

File size: 38.3 KB
Line 
1! $Id: infotrac.F90 3919 2021-05-29 08:54:49Z fhourdin $
2!
3MODULE infotrac
4
5! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
6  INTEGER, SAVE :: nqtot
7!CR: on ajoute le nombre de traceurs de l eau
8  INTEGER, SAVE :: nqo
9
10! nbtr : number of tracers not including higher order of moment or water vapor or liquid
11!        number of tracers used in the physics
12  INTEGER, SAVE :: nbtr
13
14! CRisi: on retranche les isotopes des traceurs habituels
15! On fait un tableaux d'indices des traceurs qui passeront dans phytrac
16  INTEGER, SAVE :: nqtottr
17  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: itr_indice
18
19! CRisi: nb traceurs pères= directement advectés par l'air
20  INTEGER, SAVE :: nqperes
21
22! ThL: nb traceurs INCA
23  INTEGER, SAVE :: nqINCA
24
25! ThL: nb traceurs CO2
26  INTEGER, SAVE :: nqCO2
27
28! Name variables
29  CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
30  CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics
31
32! iadv  : index of trasport schema for each tracer
33  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iadv
34
35! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the
36!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
37  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: niadv ! equivalent dyn / physique
38
39! CRisi: tableaux de fils
40  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqfils
41  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les générations
42  INTEGER, SAVE :: nqdesc_tot
43  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE    :: iqfils
44  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iqpere
45  REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi
46  PARAMETER (qperemin=1e-16,masseqmin=1e-16,ratiomin=1e-16) ! MVals
47
48! conv_flg(it)=0 : convection desactivated for tracer number it
49  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
50! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it
51  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
52
53  CHARACTER(len=4),SAVE :: type_trac
54  CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
55
56! CRisi: cas particulier des isotopes
57  LOGICAL,SAVE :: ok_isotopes,ok_iso_verif,ok_isotrac,ok_init_iso
58  INTEGER :: niso_possibles   
59  PARAMETER ( niso_possibles=5)
60  REAL, DIMENSION (niso_possibles),SAVE :: tnat,alpha_ideal
61  LOGICAL, DIMENSION(niso_possibles),SAVE ::  use_iso
62  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  iqiso ! donne indice iq en fn de (ixt,phase)
63  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_num ! donne numéro iso entre 1 et niso_possibles en fn de nqtot
64  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  iso_indnum ! donne numéro iso entre 1 et niso effectif en fn de nqtot
65  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  zone_num ! donne numéro de la zone de tracage en fn de nqtot
66  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE ::  phase_num ! donne numéro de la zone de tracage en fn de nqtot
67  INTEGER, DIMENSION(niso_possibles), SAVE :: indnum_fn_num ! donne indice entre entre 1 et niso en fonction du numéro d isotope entre 1 et niso_possibles
68  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE ::  index_trac ! numéro ixt en fn izone, indnum entre 1 et niso
69  INTEGER,SAVE :: niso,ntraceurs_zone,ntraciso
70
71#ifdef CPP_StratAer
72!--CK/OB for stratospheric aerosols
73  INTEGER, SAVE :: nbtr_bin
74  INTEGER, SAVE :: nbtr_sulgas
75  INTEGER, SAVE :: id_OCS_strat
76  INTEGER, SAVE :: id_SO2_strat
77  INTEGER, SAVE :: id_H2SO4_strat
78  INTEGER, SAVE :: id_BIN01_strat
79  INTEGER, SAVE :: id_TEST_strat
80#endif
81 
82CONTAINS
83
84  SUBROUTINE infotrac_init
85    USE control_mod, ONLY: planet_type, config_inca
86#ifdef REPROBUS
87    USE CHEM_REP, ONLY : Init_chem_rep_trac
88#endif
89    IMPLICIT NONE
90!=======================================================================
91!
92!   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
93!   -------
94!   Modif special traceur F.Forget 05/94
95!   Modif M-A Filiberti 02/02 lecture de traceur.def
96!
97!   Objet:
98!   ------
99!   GCM LMD nouvelle grille
100!
101!=======================================================================
102!   ... modification de l'integration de q ( 26/04/94 ) ....
103!-----------------------------------------------------------------------
104! Declarations
105
106    INCLUDE "dimensions.h"
107    INCLUDE "iniprint.h"
108
109! Local variables
110    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv  ! index of horizontal trasport schema
111    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv  ! index of vertical trasport schema
112
113    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca  ! index of horizontal trasport schema
114    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca  ! index of vertical trasport schema
115
116    INTEGER, ALLOCATABLE, DIMENSION(:) :: conv_flg_inca
117    INTEGER, ALLOCATABLE, DIMENSION(:) :: pbl_flg_inca
118    CHARACTER(len=8), ALLOCATABLE, DIMENSION(:) :: solsym_inca
119
120    CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
121    CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
122    CHARACTER(len=3), DIMENSION(30) :: descrq
123    CHARACTER(len=1), DIMENSION(3)  :: txts
124    CHARACTER(len=2), DIMENSION(9)  :: txtp
125    CHARACTER(len=23)               :: str1,str2
126 
127    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
128    INTEGER :: iq, new_iq, iiq, jq, ierr,itr
129    INTEGER :: ifils,ipere,generation ! CRisi
130    LOGICAL :: continu,nouveau_traceurdef
131    INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
132    CHARACTER(len=30) :: tchaine   
133
134    character(len=*),parameter :: modname="infotrac_init"
135
136!-----------------------------------------------------------------------
137! Initialization :
138!
139    txts=(/'x','y','z'/)
140    txtp=(/'x ','y ','z ','xx','xy','xz','yy','yz','zz'/)
141
142    descrq(14)='VLH'
143    descrq(10)='VL1'
144    descrq(11)='VLP'
145    descrq(12)='FH1'
146    descrq(13)='FH2'
147    descrq(16)='PPM'
148    descrq(17)='PPS'
149    descrq(18)='PPP'
150    descrq(20)='SLP'
151    descrq(30)='PRA'
152   
153
154    ! Coherence test between parameter type_trac, config_inca and preprocessing keys
155    IF (type_trac=='inca') THEN
156       WRITE(lunout,*) 'You have chosen to couple with INCA chemistry model : type_trac=', &
157            type_trac,' config_inca=',config_inca
158       IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
159          WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
160          CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
161       ENDIF
162#ifndef INCA
163       WRITE(lunout,*) 'To run this option you must add cpp key INCA and compile with INCA code'
164       CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
165#endif
166    ELSE IF (type_trac=='repr') THEN
167       WRITE(lunout,*) 'You have chosen to couple with REPROBUS chemestry model : type_trac=', type_trac
168#ifndef REPROBUS
169       WRITE(lunout,*) 'To run this option you must add cpp key REPROBUS and compile with REPRPBUS code'
170       CALL abort_gcm('infotrac_init','You must compile with cpp key REPROBUS',1)
171#endif
172    ELSE IF (type_trac == 'co2i') THEN
173       WRITE(lunout,*) 'You have chosen to run with CO2 cycle: type_trac=', type_trac
174    ELSE IF (type_trac == 'coag') THEN
175       WRITE(lunout,*) 'Tracers are treated for COAGULATION tests : type_trac=', type_trac
176#ifndef CPP_StratAer
177       WRITE(lunout,*) 'To run this option you must add cpp key StratAer and compile with StratAer code'
178       CALL abort_gcm('infotrac_init','You must compile with cpp key StratAer',1)
179#endif
180    ELSE IF (type_trac == 'lmdz') THEN
181       WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac
182    ELSE IF (type_trac == 'inco') THEN ! ThL
183       WRITE(lunout,*) 'Using jointly INCA and CO2 cycle: type_trac =', type_trac
184       IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
185          WRITE(lunout,*) 'Incoherence between type_trac and config_inca. Model stops. Modify run.def'
186          CALL abort_gcm('infotrac_init','Incoherence between type_trac and config_inca',1)
187       ENDIF
188#ifndef INCA
189       WRITE(lunout,*) 'To run this option you must add cpp key INCA and compilewith INCA code'
190       CALL abort_gcm('infotrac_init','You must compile with cpp key INCA',1)
191#endif   
192    ELSE
193       WRITE(lunout,*) 'type_trac=',type_trac,' not possible. Model stops'
194       CALL abort_gcm('infotrac_init','bad parameter',1)
195    ENDIF
196
197    ! Test if config_inca is other then none for run without INCA
198    IF (type_trac/='inca' .AND. type_trac/='inco' .AND. config_inca/='none') THEN
199       WRITE(lunout,*) 'config_inca will now be changed to none as you do not couple with INCA model'
200       config_inca='none'
201    ENDIF
202
203!-----------------------------------------------------------------------
204!
205! 1) Get the true number of tracers + water vapor/liquid
206!    Here true tracers (nqtrue) means declared tracers (only first order)
207!
208!-----------------------------------------------------------------------
209    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
210       IF (type_trac=='co2i') THEN
211          nqCO2 = 1
212       ELSE
213          nqCO2 = 0
214       ENDIF
215       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
216       IF(ierr.EQ.0) THEN
217          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
218          READ(90,*) nqtrue
219          write(lunout,*) 'nqtrue=',nqtrue
220       ELSE
221          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
222          WRITE(lunout,*) trim(modname),': WARNING using defaut values'
223          IF (planet_type=='earth') THEN
224            nqtrue=4 ! Default value for Earth
225          ELSE
226            nqtrue=1 ! Default value for other planets
227          ENDIF
228       ENDIF
229!jyg<
230!!       if ( planet_type=='earth') then
231!!         ! For Earth, water vapour & liquid tracers are not in the physics
232!!         nbtr=nqtrue-2
233!!       else
234!!         ! Other planets (for now); we have the same number of tracers
235!!         ! in the dynamics than in the physics
236!!         nbtr=nqtrue
237!!       endif
238!>jyg
239    ELSE ! type_trac=inca or inco
240       IF (type_trac=='inco') THEN
241          nqCO2 = 1
242       ELSE
243          nqCO2 = 0
244       ENDIF
245!jyg<
246       ! The traceur.def file is used to define the number "nqo" of water phases
247       ! present in the simulation. Default : nqo = 2.
248       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
249       IF(ierr.EQ.0) THEN
250          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
251          READ(90,*) nqo
252       ELSE
253          WRITE(lunout,*) trim(modname),': Using default value for nqo'
254          nqo=2
255       ENDIF
256       IF (nqo /= 2 .AND. nqo /= 3 ) THEN
257          IF (nqo == 4 .AND. type_trac=='inco') THEN ! ThL
258             WRITE(lunout,*) trim(modname),': you are coupling with INCA, and also using CO2i.'
259             nqo = 3    ! A améliorier... je force 3 traceurs eau...  ThL
260             WRITE(lunout,*) trim(modname),': nqo = ',nqo
261          ELSE
262          WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowed. Only 2 or 3 water phases allowed'
263          CALL abort_gcm('infotrac_init','Bad number of water phases',1)
264          ENDIF
265       ENDIF
266       ! nbtr has been read from INCA by init_const_lmdz() in gcm.F
267#ifdef INCA
268       CALL Init_chem_inca_trac(nqINCA)
269#else
270       nqINCA=0
271#endif
272       nbtr=nqINCA+nqCO2
273       nqtrue=nbtr+nqo
274       WRITE(lunout,*) trim(modname),': nqo = ',nqo
275       WRITE(lunout,*) trim(modname),': nbtr = ',nbtr
276       WRITE(lunout,*) trim(modname),': nqtrue = ',nqtrue
277       WRITE(lunout,*) trim(modname),': nqCO2 = ',nqCO2
278       WRITE(lunout,*) trim(modname),': nqINCA = ',nqINCA
279       ALLOCATE(hadv_inca(nqINCA), vadv_inca(nqINCA), conv_flg_inca(nqINCA), pbl_flg_inca(nqINCA), solsym_inca(nqINCA))
280    ENDIF   ! type_trac 'inca' ou 'inco'
281!>jyg
282
283    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
284       WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowed. 2 tracers is the minimum'
285       CALL abort_gcm('infotrac_init','Not enough tracers',1)
286    ENDIF
287   
288!jyg<
289! Transfert number of tracers to Reprobus
290!!    IF (type_trac == 'repr') THEN
291!!#ifdef REPROBUS
292!!       CALL Init_chem_rep_trac(nbtr)
293!!#endif
294!!    ENDIF
295!>jyg
296       
297!
298! Allocate variables depending on nqtrue
299!
300    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
301
302!
303!jyg<
304!!    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
305!!    conv_flg(:) = 1 ! convection activated for all tracers
306!!    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
307!>jyg
308
309!-----------------------------------------------------------------------
310! 2)     Choix  des schemas d'advection pour l'eau et les traceurs
311!
312!     iadv = 1    schema  transport type "humidite specifique LMD"
313!     iadv = 2    schema   amont
314!     iadv = 14   schema  Van-leer + humidite specifique
315!                            Modif F.Codron
316!     iadv = 10   schema  Van-leer (retenu pour l'eau vapeur et liquide)
317!     iadv = 11   schema  Van-Leer pour hadv et version PPM (Monotone) pour vadv
318!     iadv = 12   schema  Frederic Hourdin I
319!     iadv = 13   schema  Frederic Hourdin II
320!     iadv = 16   schema  PPM Monotone(Collela & Woodward 1984)
321!     iadv = 17   schema  PPM Semi Monotone (overshoots autorisés)
322!     iadv = 18   schema  PPM Positif Defini (overshoots undershoots autorisés)
323!     iadv = 20   schema  Slopes
324!     iadv = 30   schema  Prather
325!
326!        Dans le tableau q(ij,l,iq) : iq = 1  pour l'eau vapeur
327!                                     iq = 2  pour l'eau liquide
328!       Et eventuellement             iq = 3,nqtot pour les autres traceurs
329!
330!        iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
331!------------------------------------------------------------------------
332!
333!    Get choice of advection schema from file tracer.def or from INCA
334!---------------------------------------------------------------------
335    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
336       IF(ierr.EQ.0) THEN
337          ! Continue to read tracer.def
338          DO iq=1,nqtrue
339
340             write(*,*) 'infotrac 237: iq=',iq
341             ! CRisi: ajout du nom du fluide transporteur
342             ! mais rester retro compatible
343             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
344             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
345             write(lunout,*) 'tchaine=',trim(tchaine)
346             write(*,*) 'infotrac 238: IOstatus=',IOstatus
347             if (IOstatus.ne.0) then
348                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
349             endif
350             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
351             ! espace ou pas au milieu de la chaine.
352             continu=.true.
353             nouveau_traceurdef=.false.
354             iiq=1
355             do while (continu)
356                if (tchaine(iiq:iiq).eq.' ') then
357                  nouveau_traceurdef=.true.
358                  continu=.false.
359                else if (iiq.lt.LEN_TRIM(tchaine)) then
360                  iiq=iiq+1
361                else
362                  continu=.false.
363                endif
364             enddo
365             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
366             if (nouveau_traceurdef) then
367                write(lunout,*) 'C''est la nouvelle version de traceur.def'
368                tnom_0(iq)=tchaine(1:iiq-1)
369                tnom_transp(iq)=tchaine(iiq+1:30)
370             else
371                write(lunout,*) 'C''est l''ancienne version de traceur.def'
372                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
373                tnom_0(iq)=tchaine
374                tnom_transp(iq) = 'air'
375             endif
376             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
377             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
378
379          ENDDO!DO iq=1,nqtrue
380          CLOSE(90) 
381
382       ELSE ! Without tracer.def, set default values
383         if (planet_type=="earth") then
384          ! for Earth, default is to have 4 tracers
385          hadv(1) = 14
386          vadv(1) = 14
387          tnom_0(1) = 'H2Ov'
388          tnom_transp(1) = 'air'
389          hadv(2) = 10
390          vadv(2) = 10
391          tnom_0(2) = 'H2Ol'
392          tnom_transp(2) = 'air'
393          hadv(3) = 10
394          vadv(3) = 10
395          tnom_0(3) = 'RN'
396          tnom_transp(3) = 'air'
397          hadv(4) = 10
398          vadv(4) = 10
399          tnom_0(4) = 'PB'
400          tnom_transp(4) = 'air'
401         else ! default for other planets
402          hadv(1) = 10
403          vadv(1) = 10
404          tnom_0(1) = 'dummy'
405          tnom_transp(1) = 'dummy'
406         endif ! of if (planet_type=="earth")
407       ENDIF
408       
409       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
410       WRITE(lunout,*) trim(modname),': nombre total de traceurs ',nqtrue
411       DO iq=1,nqtrue
412          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
413       END DO
414
415       IF ( planet_type=='earth') THEN
416         !CR: nombre de traceurs de l eau
417         IF (tnom_0(3) == 'H2Oi') THEN
418            nqo=3
419         ELSE
420            nqo=2
421         ENDIF
422         ! For Earth, water vapour & liquid tracers are not in the physics
423         nbtr=nqtrue-nqo
424       ELSE
425         ! Other planets (for now); we have the same number of tracers
426         ! in the dynamics than in the physics
427         nbtr=nqtrue
428       ENDIF
429
430#ifdef CPP_StratAer
431       IF (type_trac == 'coag') THEN
432         nbtr_bin=0
433         nbtr_sulgas=0
434         DO iq=1,nqtrue
435           IF (tnom_0(iq)(1:3)=='BIN') THEN !check if tracer name contains 'BIN'
436             nbtr_bin=nbtr_bin+1
437           ENDIF
438           IF (tnom_0(iq)(1:3)=='GAS') THEN !check if tracer name contains 'GAS'
439             nbtr_sulgas=nbtr_sulgas+1
440           ENDIF
441         ENDDO
442         print*,'nbtr_bin=',nbtr_bin
443         print*,'nbtr_sulgas=',nbtr_sulgas
444         DO iq=1,nqtrue
445           IF (tnom_0(iq)=='GASOCS') THEN
446             id_OCS_strat=iq-nqo
447           ENDIF
448           IF (tnom_0(iq)=='GASSO2') THEN
449             id_SO2_strat=iq-nqo
450           ENDIF
451           IF (tnom_0(iq)=='GASH2SO4') THEN
452             id_H2SO4_strat=iq-nqo
453           ENDIF
454           IF (tnom_0(iq)=='BIN01') THEN
455             id_BIN01_strat=iq-nqo
456           ENDIF
457           IF (tnom_0(iq)=='GASTEST') THEN
458             id_TEST_strat=iq-nqo
459           ENDIF
460         ENDDO
461         print*,'id_OCS_strat  =',id_OCS_strat
462         print*,'id_SO2_strat  =',id_SO2_strat
463         print*,'id_H2SO4_strat=',id_H2SO4_strat
464         print*,'id_BIN01_strat=',id_BIN01_strat
465       ENDIF
466#endif
467
468    ENDIF  ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag' .OR. type_trac = 'co2i')
469!jyg<
470!
471
472!-----------------------------------------------------------------------
473!
474! 3) Verify if advection schema 20 or 30 choosen
475!    Calculate total number of tracers needed: nqtot
476!    Allocate variables depending on total number of tracers
477!-----------------------------------------------------------------------
478    new_iq=0
479    DO iq=1,nqtrue
480       ! Add tracers for certain advection schema
481       IF (hadv(iq)<20 .AND. vadv(iq)<20 ) THEN
482          new_iq=new_iq+1  ! no tracers added
483       ELSE IF (hadv(iq)==20 .AND. vadv(iq)==20 ) THEN
484          new_iq=new_iq+4  ! 3 tracers added
485       ELSE IF (hadv(iq)==30 .AND. vadv(iq)==30 ) THEN
486          new_iq=new_iq+10 ! 9 tracers added
487       ELSE
488          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
489          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
490       ENDIF
491    END DO
492   
493    IF (new_iq /= nqtrue) THEN
494       ! The choice of advection schema imposes more tracers
495       ! Assigne total number of tracers
496       nqtot = new_iq
497
498       WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
499       WRITE(lunout,*) 'makes it necessary to add tracers'
500       WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers'
501       WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed'
502
503    ELSE
504       ! The true number of tracers is also the total number
505       nqtot = nqtrue
506    ENDIF
507
508!
509! Allocate variables with total number of tracers, nqtot
510!
511    ALLOCATE(tname(nqtot), ttext(nqtot))
512    ALLOCATE(iadv(nqtot), niadv(nqtot))
513
514!-----------------------------------------------------------------------
515!
516! 4) Determine iadv, long and short name
517!
518!-----------------------------------------------------------------------
519    new_iq=0
520    DO iq=1,nqtrue
521       new_iq=new_iq+1
522
523       ! Verify choice of advection schema
524       IF (hadv(iq)==vadv(iq)) THEN
525          iadv(new_iq)=hadv(iq)
526       ELSE IF (hadv(iq)==10 .AND. vadv(iq)==16) THEN
527          iadv(new_iq)=11
528       ELSE
529          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
530
531          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
532       ENDIF
533     
534       str1=tnom_0(iq)
535       tname(new_iq)= tnom_0(iq)
536       IF (iadv(new_iq)==0) THEN
537          ttext(new_iq)=trim(str1)
538       ELSE
539          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
540       ENDIF
541
542       ! schemas tenant compte des moments d'ordre superieur
543       str2=ttext(new_iq)
544       IF (iadv(new_iq)==20) THEN
545          DO jq=1,3
546             new_iq=new_iq+1
547             iadv(new_iq)=-20
548             ttext(new_iq)=trim(str2)//txts(jq)
549             tname(new_iq)=trim(str1)//txts(jq)
550          END DO
551       ELSE IF (iadv(new_iq)==30) THEN
552          DO jq=1,9
553             new_iq=new_iq+1
554             iadv(new_iq)=-30
555             ttext(new_iq)=trim(str2)//txtp(jq)
556             tname(new_iq)=trim(str1)//txtp(jq)
557          END DO
558       ENDIF
559    END DO
560
561!
562! Find vector keeping the correspodence between true and total tracers
563!
564    niadv(:)=0
565    iiq=0
566    DO iq=1,nqtot
567       IF(iadv(iq).GE.0) THEN
568          ! True tracer
569          iiq=iiq+1
570          niadv(iiq)=iq
571       ENDIF
572    END DO
573
574
575    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
576    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
577    DO iq=1,nqtot
578       WRITE(lunout,*) iadv(iq),niadv(iq),&
579       ' ',trim(tname(iq)),' ',trim(ttext(iq))
580    END DO
581
582!
583! Test for advection schema.
584! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) .
585!
586    DO iq=1,nqtot
587       IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
588          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
589          CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
590       ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
591          WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
592          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
593       ENDIF
594    END DO
595
596
597! CRisi: quels sont les traceurs fils et les traceurs pères.
598! initialiser tous les tableaux d'indices liés aux traceurs familiaux
599! + vérifier que tous les pères sont écrits en premières positions
600    ALLOCATE(nqfils(nqtot),nqdesc(nqtot))   
601    ALLOCATE(iqfils(nqtot,nqtot))   
602    ALLOCATE(iqpere(nqtot))
603    nqperes=0
604    nqfils(:)=0
605    nqdesc(:)=0
606    iqfils(:,:)=0
607    iqpere(:)=0
608    nqdesc_tot=0   
609    DO iq=1,nqtot
610      if (tnom_transp(iq) == 'air') then
611        ! ceci est un traceur père
612        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
613        nqperes=nqperes+1
614        iqpere(iq)=0
615      else !if (tnom_transp(iq) == 'air') then
616        ! ceci est un fils. Qui est son père?
617        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils'
618        continu=.true.
619        ipere=1
620        do while (continu)           
621          if (tnom_transp(iq) == tnom_0(ipere)) then
622            ! Son père est ipere
623            WRITE(lunout,*) 'Le traceur',iq,'appele ', &
624      &          trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
625            if (iq.eq.ipere) then
626                CALL abort_gcm('infotrac_init','Un fils est son propre pere',1)
627            endif
628            nqfils(ipere)=nqfils(ipere)+1 
629            iqfils(nqfils(ipere),ipere)=iq
630            iqpere(iq)=ipere         
631            continu=.false.
632          else !if (tnom_transp(iq) == tnom_0(ipere)) then
633            ipere=ipere+1
634            if (ipere.gt.nqtot) then
635                WRITE(lunout,*) 'Le traceur',iq,'appele ', &
636      &          trim(tnom_0(iq)),', est orphelin.'
637                CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
638            endif !if (ipere.gt.nqtot) then
639          endif !if (tnom_transp(iq) == tnom_0(ipere)) then
640        enddo !do while (continu)
641      endif !if (tnom_transp(iq) == 'air') then
642    enddo !DO iq=1,nqtot
643    WRITE(lunout,*) 'infotrac: nqperes=',nqperes   
644    WRITE(lunout,*) 'nqfils=',nqfils
645    WRITE(lunout,*) 'iqpere=',iqpere
646    WRITE(lunout,*) 'iqfils=',iqfils
647
648! Calculer le nombre de descendants à partir de iqfils et de nbfils
649    DO iq=1,nqtot   
650      generation=0
651      continu=.true.
652      ifils=iq
653      do while (continu)
654        ipere=iqpere(ifils)
655        if (ipere.gt.0) then
656         nqdesc(ipere)=nqdesc(ipere)+1   
657         nqdesc_tot=nqdesc_tot+1     
658         iqfils(nqdesc(ipere),ipere)=iq
659         ifils=ipere
660         generation=generation+1
661        else !if (ipere.gt.0) then
662         continu=.false.
663        endif !if (ipere.gt.0) then
664      enddo !do while (continu)   
665      WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation
666    enddo !DO iq=1,nqtot
667    WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc
668    WRITE(lunout,*) 'iqfils=',iqfils
669    WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot
670
671! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas
672! que 10 et 14 si des pères ont des fils
673    do iq=1,nqtot
674      if (iqpere(iq).gt.0) then
675        ! ce traceur a un père qui n'est pas l'air
676        ! Seul le schéma 10 est autorisé
677        if (iadv(iq)/=10) then
678           WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons'
679          CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1)
680        endif
681        ! Le traceur père ne peut être advecté que par schéma 10 ou 14:
682        IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
683          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers'
684          CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1)
685        endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
686     endif !if (iqpere(iq).gt.0) the
687    enddo !do iq=1,nqtot
688
689
690
691! detecter quels sont les traceurs isotopiques parmi des traceurs
692    call infotrac_isoinit(tnom_0,nqtrue)
693
694!    if (ntraciso.gt.0) then
695! le 18 sep 2020: on enlève la condition ntraciso.gt.0 car nqtottr doit être
696! connu même si il n'y a pas d'isotopes!
697        write(lunout,*) 'infotrac 702: nbtr,ntraciso=',nbtr,ntraciso
698! retrancher les traceurs isotopiques de la liste des traceurs qui passent dans
699! phytrac
700        nbtr=nbtr-nqo*ntraciso
701
702! faire un tableau d'indice des traceurs qui passeront dans phytrac
703        nqtottr=nqtot-nqo*(1+ntraciso)
704        write(lunout,*) 'infotrac 704: nqtottr,nqtot,nqo=',nqtottr,nqtot,nqo
705        ! Rq: nqtottr n'est pas forcément égal à nbtr dans le cas où new_iq /= nqtrue
706        ALLOCATE (itr_indice(nqtottr)) 
707        itr_indice(:)=0 
708        itr=0
709        do iq=nqo+1, nqtot
710          if (iso_num(iq).eq.0) then
711            itr=itr+1
712            write(*,*) 'itr=',itr
713            itr_indice(itr)=iq
714          endif !if (iso_num(iq).eq.0) then
715        enddo
716        if (itr.ne.nqtottr) then
717            CALL abort_gcm('infotrac_init','pb dans le calcul de nqtottr',1)
718        endif
719        write(lunout,*) 'itr_indice=',itr_indice
720!    endif !if (ntraciso.gt.0) then
721
722
723! Transfert number of tracers to Reprobus
724    IF (type_trac == 'repr') THEN
725#ifdef REPROBUS
726       CALL Init_chem_rep_trac(nbtr,nqo,tnom_0)
727#endif
728    ENDIF
729!
730! Allocate variables depending on nbtr
731!
732    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
733    conv_flg(:) = 1 ! convection activated for all tracers
734    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
735!
736!!    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
737!
738    IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN   ! config_inca='aero' ou 'chem'
739!>jyg
740! le module de chimie fournit les noms des traceurs
741! et les schemas d'advection associes. excepte pour ceux lus
742! dans traceur.def
743       IF (ierr .eq. 0) then
744          DO iq=1,nqo+nqCO2
745
746             write(*,*) 'infotrac 237: iq=',iq
747             ! CRisi: ajout du nom du fluide transporteur
748             ! mais rester retro compatible
749             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
750             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
751             write(lunout,*) 'tchaine=',trim(tchaine)
752             write(*,*) 'infotrac 238: IOstatus=',IOstatus
753             if (IOstatus.ne.0) then
754                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
755             endif
756             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
757             ! espace ou pas au milieu de la chaine.
758             continu=.true.
759             nouveau_traceurdef=.false.
760             iiq=1
761             do while (continu)
762                if (tchaine(iiq:iiq).eq.' ') then
763                  nouveau_traceurdef=.true.
764                  continu=.false.
765                else if (iiq.lt.LEN_TRIM(tchaine)) then
766                  iiq=iiq+1
767                else
768                  continu=.false.
769                endif
770             enddo
771             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
772             if (nouveau_traceurdef) then
773                write(lunout,*) 'C''est la nouvelle version de traceur.def'
774                tnom_0(iq)=tchaine(1:iiq-1)
775                tnom_transp(iq)=tchaine(iiq+1:30)
776             else
777                write(lunout,*) 'C''est l''ancienne version de traceur.def'
778                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
779                tnom_0(iq)=tchaine
780                tnom_transp(iq) = 'air'
781             endif
782             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
783             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
784
785          ENDDO  !DO iq=1,nqo
786          CLOSE(90) 
787       ELSE  !! if traceur.def doesn't exist
788          tnom_0(1)='H2Ov'
789          tnom_transp(1) = 'air'
790          tnom_0(2)='H2Ol'
791          tnom_transp(2) = 'air'
792          hadv(1) = 10
793          hadv(2) = 10
794          vadv(1) = 10
795          vadv(2) = 10
796       ENDIF
797 
798#ifdef INCA
799       CALL init_transport( &
800            hadv_inca, &
801            vadv_inca, &
802            conv_flg_inca, &
803            pbl_flg_inca,  &
804            solsym_inca)
805
806       conv_flg(1+nqCO2:nbtr) = conv_flg_inca
807       pbl_flg(1+nqCO2:nbtr) = pbl_flg_inca
808       solsym(1+nqCO2:nbtr) = solsym_inca
809
810       IF (type_trac == 'inco') THEN
811          conv_flg(1:nqCO2) = 1
812          pbl_flg(1:nqCO2) = 1
813          solsym(1:nqCO2) = 'CO2'
814       ENDIF
815#endif
816
817!jyg<
818       DO iq = nqo+nqCO2+1, nqtrue
819          hadv(iq) = hadv_inca(iq-nqo-nqCO2)
820          vadv(iq) = vadv_inca(iq-nqo-nqCO2)
821          tnom_0(iq)=solsym_inca(iq-nqo-nqCO2)
822          tnom_transp(iq) = 'air'
823       END DO
824
825    ENDIF ! (type_trac == 'inca' or 'inco')
826
827!-----------------------------------------------------------------------
828! Finalize :
829!
830    DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
831
832    WRITE(lunout,*) 'infotrac init fin'
833
834  END SUBROUTINE infotrac_init
835
836  SUBROUTINE infotrac_isoinit(tnom_0,nqtrue)
837
838#ifdef CPP_IOIPSL
839  use IOIPSL
840#else
841  ! if not using IOIPSL, we still need to use (a local version of) getin
842  use ioipsl_getincom
843#endif
844  implicit none
845 
846    ! inputs
847    INTEGER nqtrue
848    CHARACTER(len=30) tnom_0(nqtrue)
849   
850    ! locals   
851    CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso
852    INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso
853    INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind
854    INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone
855    CHARACTER(len=19) :: tnom_trac
856    INCLUDE "iniprint.h"
857
858    tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
859
860    ALLOCATE(nb_iso(niso_possibles,nqo))
861    ALLOCATE(nb_isoind(nqo))
862    ALLOCATE(nb_traciso(niso_possibles,nqo))
863    ALLOCATE(iso_num(nqtot))
864    ALLOCATE(iso_indnum(nqtot))
865    ALLOCATE(zone_num(nqtot))
866    ALLOCATE(phase_num(nqtot))
867     
868    iso_num(:)=0
869    iso_indnum(:)=0
870    zone_num(:)=0
871    phase_num(:)=0
872    indnum_fn_num(:)=0
873    use_iso(:)=.false. 
874    nb_iso(:,:)=0 
875    nb_isoind(:)=0     
876    nb_traciso(:,:)=0
877    niso=0
878    ntraceurs_zone=0 
879    ntraceurs_zone_prec=0
880    ntraciso=0
881
882    do iq=nqo+1,nqtot
883!       write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
884       do phase=1,nqo   
885        do ixt= 1,niso_possibles   
886         tnom_trac=trim(tnom_0(phase))//'_'
887         tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt))
888!         write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac     
889         IF (tnom_0(iq) == tnom_trac) then
890!          write(lunout,*) 'Ce traceur est un isotope'
891          nb_iso(ixt,phase)=nb_iso(ixt,phase)+1   
892          nb_isoind(phase)=nb_isoind(phase)+1   
893          iso_num(iq)=ixt
894          iso_indnum(iq)=nb_isoind(phase)
895          indnum_fn_num(ixt)=iso_indnum(iq)
896          phase_num(iq)=phase
897!          write(lunout,*) 'iso_num(iq)=',iso_num(iq)
898!          write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq)
899!          write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt)
900!          write(lunout,*) 'phase_num(iq)=',phase_num(iq)
901          goto 20
902         else if (iqpere(iq).gt.0) then         
903          if (tnom_0(iqpere(iq)) == tnom_trac) then
904!           write(lunout,*) 'Ce traceur est le fils d''un isotope'
905           ! c'est un traceur d'isotope
906           nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1
907           iso_num(iq)=ixt
908           iso_indnum(iq)=indnum_fn_num(ixt)
909           zone_num(iq)=nb_traciso(ixt,phase)
910           phase_num(iq)=phase
911!           write(lunout,*) 'iso_num(iq)=',iso_num(iq)
912!           write(lunout,*) 'phase_num(iq)=',phase_num(iq)
913!           write(lunout,*) 'zone_num(iq)=',zone_num(iq)
914           goto 20
915          endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
916         endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
917        enddo !do ixt= niso_possibles
918       enddo !do phase=1,nqo
919  20   continue
920      enddo !do iq=1,nqtot
921
922!      write(lunout,*) 'iso_num=',iso_num
923!      write(lunout,*) 'iso_indnum=',iso_indnum
924!      write(lunout,*) 'zone_num=',zone_num 
925!      write(lunout,*) 'phase_num=',phase_num
926!      write(lunout,*) 'indnum_fn_num=',indnum_fn_num
927
928      do ixt= 1,niso_possibles 
929
930        if (nb_iso(ixt,1).eq.1) then
931          ! on vérifie que toutes les phases ont le même nombre de
932          ! traceurs
933          do phase=2,nqo
934            if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
935!              write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
936              CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1)
937            endif
938          enddo !do phase=2,nqo
939
940          niso=niso+1
941          use_iso(ixt)=.true.
942          ntraceurs_zone=nb_traciso(ixt,1)
943
944          ! on vérifie que toutes les phases ont le même nombre de
945          ! traceurs
946          do phase=2,nqo
947            if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then
948              write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase)
949              write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone
950              CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
951            endif 
952          enddo  !do phase=2,nqo
953          ! on vérifie que tous les isotopes ont le même nombre de
954          ! traceurs
955          if (ntraceurs_zone_prec.gt.0) then               
956            if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
957              ntraceurs_zone_prec=ntraceurs_zone
958            else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
959              write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone   
960              CALL abort_gcm('infotrac_init', &
961               &'Isotope tracers are not well defined in traceur.def',1)           
962            endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
963           endif !if (ntraceurs_zone_prec.gt.0) then
964
965        else if (nb_iso(ixt,1).ne.0) then
966           WRITE(lunout,*) 'nqo,ixt=',nqo,ixt
967           WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1)   
968           CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1)     
969        endif   !if (nb_iso(ixt,1).eq.1) then       
970    enddo ! do ixt= niso_possibles
971
972    ! dimensions isotopique:
973    ntraciso=niso*(ntraceurs_zone+1)
974!    WRITE(lunout,*) 'niso=',niso
975!    WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso   
976 
977    ! flags isotopiques:
978    if (niso.gt.0) then
979        ok_isotopes=.true.
980    else
981        ok_isotopes=.false.
982    endif
983!    WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
984 
985    if (ok_isotopes) then
986        ok_iso_verif=.false.
987        call getin('ok_iso_verif',ok_iso_verif)
988        ok_init_iso=.false.
989        call getin('ok_init_iso',ok_init_iso)
990        tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/)
991        alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/)
992    endif !if (ok_isotopes) then 
993!    WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
994!    WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
995
996    if (ntraceurs_zone.gt.0) then
997        ok_isotrac=.true.
998    else
999        ok_isotrac=.false.
1000    endif   
1001!    WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
1002
1003    ! remplissage du tableau iqiso(ntraciso,phase)
1004    ALLOCATE(iqiso(ntraciso,nqo))   
1005    iqiso(:,:)=0     
1006    do iq=1,nqtot
1007        if (iso_num(iq).gt.0) then
1008          ixt=iso_indnum(iq)+zone_num(iq)*niso
1009          iqiso(ixt,phase_num(iq))=iq
1010        endif
1011    enddo
1012!    WRITE(lunout,*) 'iqiso=',iqiso
1013
1014    ! replissage du tableau index_trac(ntraceurs_zone,niso)
1015    ALLOCATE(index_trac(ntraceurs_zone,niso)) 
1016    if (ok_isotrac) then
1017        do iiso=1,niso
1018          do izone=1,ntraceurs_zone
1019             index_trac(izone,iiso)=iiso+izone*niso
1020          enddo
1021        enddo
1022    else !if (ok_isotrac) then     
1023        index_trac(:,:)=0.0
1024    endif !if (ok_isotrac) then
1025!    write(lunout,*) 'index_trac=',index_trac   
1026
1027! Finalize :
1028    DEALLOCATE(nb_iso)
1029
1030  END SUBROUTINE infotrac_isoinit
1031
1032END MODULE infotrac
Note: See TracBrowser for help on using the repository browser.