source: LMDZ6/trunk/libf/dyn3d_common/infotrac.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by lmdz-users, 3 years ago

Modifications from Thibaut to create an ESM with interactive CO2 + INCA aerosols

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