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

Last change on this file since 3923 was 3923, checked in by fhourdin, 3 years ago

Version modifiee par Camille pour compatibilité avec la physique isotopique.

  • 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: 38.3 KB
RevLine 
[1279]1! $Id: infotrac.F90 3923 2021-06-08 09:31:06Z fhourdin $
2!
[1114]3MODULE infotrac
4
5! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
6  INTEGER, SAVE :: nqtot
[2086]7!CR: on ajoute le nombre de traceurs de l eau
8  INTEGER, SAVE :: nqo
[1114]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
[3923]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
[2270]20  INTEGER, SAVE :: nqperes
21
[3870]22! ThL: nb traceurs INCA
[3865]23  INTEGER, SAVE :: nqINCA
24
[3870]25! ThL: nb traceurs CO2
26  INTEGER, SAVE :: nqCO2
27
[1114]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
[2270]39! CRisi: tableaux de fils
40  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqfils
[3923]41  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: nqdesc ! nombres de fils + nombre de tous les petits fils sur toutes les g?n?rations
[2270]42  INTEGER, SAVE :: nqdesc_tot
43  INTEGER, ALLOCATABLE, DIMENSION(:,:), SAVE    :: iqfils
44  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE    :: iqpere
[3800]45  REAL :: qperemin,masseqmin,ratiomin ! MVals et CRisi
46  PARAMETER (qperemin=1e-16,masseqmin=1e-16,ratiomin=1e-16) ! MVals
[2270]47
[1279]48! conv_flg(it)=0 : convection desactivated for tracer number it
[1114]49  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: conv_flg
[1279]50! pbl_flg(it)=0  : boundary layer diffusion desactivaded for tracer number it
[1114]51  INTEGER, ALLOCATABLE, DIMENSION(:), SAVE  :: pbl_flg
52
[1279]53  CHARACTER(len=4),SAVE :: type_trac
[2180]54  CHARACTER(len=8),DIMENSION(:),ALLOCATABLE, SAVE :: solsym
[3870]55
[2690]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)
[3923]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
[2690]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
[1279]81 
[1114]82CONTAINS
83
84  SUBROUTINE infotrac_init
[2391]85    USE control_mod, ONLY: planet_type, config_inca
[1565]86#ifdef REPROBUS
87    USE CHEM_REP, ONLY : Init_chem_rep_trac
88#endif
[1114]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
[2566]113    INTEGER, ALLOCATABLE, DIMENSION(:) :: hadv_inca  ! index of horizontal trasport schema
114    INTEGER, ALLOCATABLE, DIMENSION(:) :: vadv_inca  ! index of vertical trasport schema
115
[3870]116    INTEGER, ALLOCATABLE, DIMENSION(:) :: conv_flg_inca
117    INTEGER, ALLOCATABLE, DIMENSION(:) :: pbl_flg_inca
118    CHARACTER(len=8), ALLOCATABLE, DIMENSION(:) :: solsym_inca
119
[3923]120    CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
121    CHARACTER(len=30), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
[1114]122    CHARACTER(len=3), DIMENSION(30) :: descrq
123    CHARACTER(len=1), DIMENSION(3)  :: txts
124    CHARACTER(len=2), DIMENSION(9)  :: txtp
[1454]125    CHARACTER(len=23)               :: str1,str2
[1114]126 
127    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
[3923]128    INTEGER :: iq, new_iq, iiq, jq, ierr,itr
[2270]129    INTEGER :: ifils,ipere,generation ! CRisi
130    LOGICAL :: continu,nouveau_traceurdef
131    INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
[3923]132    CHARACTER(len=30) :: tchaine   
[1454]133
134    character(len=*),parameter :: modname="infotrac_init"
[3865]135
[1114]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'
[1279]152   
[1114]153
[1569]154    ! Coherence test between parameter type_trac, config_inca and preprocessing keys
[1563]155    IF (type_trac=='inca') THEN
[3865]156       WRITE(lunout,*) 'You have chosen to couple with INCA chemistry model : type_trac=', &
[1563]157            type_trac,' config_inca=',config_inca
[2180]158       IF (config_inca/='aero' .AND. config_inca/='aeNP' .AND. config_inca/='chem') THEN
[1563]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)
[3872]161       ENDIF
[1569]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
[1565]166    ELSE IF (type_trac=='repr') THEN
[3865]167       WRITE(lunout,*) 'You have chosen to couple with REPROBUS chemestry model : type_trac=', type_trac
[1569]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
[3361]172    ELSE IF (type_trac == 'co2i') THEN
173       WRITE(lunout,*) 'You have chosen to run with CO2 cycle: type_trac=', type_trac
[2690]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
[1563]180    ELSE IF (type_trac == 'lmdz') THEN
181       WRITE(lunout,*) 'Tracers are treated in LMDZ only : type_trac=', type_trac
[3865]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)
[3872]187       ENDIF
[3865]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   
[1279]192    ELSE
[1563]193       WRITE(lunout,*) 'type_trac=',type_trac,' not possible. Model stops'
194       CALL abort_gcm('infotrac_init','bad parameter',1)
[3872]195    ENDIF
[1279]196
[1563]197    ! Test if config_inca is other then none for run without INCA
[3865]198    IF (type_trac/='inca' .AND. type_trac/='inco' .AND. config_inca/='none') THEN
[1563]199       WRITE(lunout,*) 'config_inca will now be changed to none as you do not couple with INCA model'
200       config_inca='none'
[3872]201    ENDIF
[1563]202
[1114]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!-----------------------------------------------------------------------
[3361]209    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
[3872]210       IF (type_trac=='co2i') THEN
[3870]211          nqCO2 = 1
212       ELSE
213          nqCO2 = 0
214       ENDIF
[1114]215       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
216       IF(ierr.EQ.0) THEN
[1454]217          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
[1114]218          READ(90,*) nqtrue
[2270]219          write(lunout,*) 'nqtrue=',nqtrue
[1114]220       ELSE
[1454]221          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
222          WRITE(lunout,*) trim(modname),': WARNING using defaut values'
[2690]223          IF (planet_type=='earth') THEN
[1454]224            nqtrue=4 ! Default value for Earth
[2690]225          ELSE
[1454]226            nqtrue=1 ! Default value for other planets
[2690]227          ENDIF
228       ENDIF
[2262]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
[3870]239    ELSE ! type_trac=inca or inco
[3872]240       IF (type_trac=='inco') THEN
[3870]241          nqCO2 = 1
242       ELSE
243          nqCO2 = 0
244       ENDIF
[2262]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
[2362]256       IF (nqo /= 2 .AND. nqo /= 3 ) THEN
[3865]257          IF (nqo == 4 .AND. type_trac=='inco') THEN ! ThL
258             WRITE(lunout,*) trim(modname),': you are coupling with INCA, and also using CO2i.'
[3923]259             nqo = 3    ! A am?liorier... je force 3 traceurs eau...  ThL
[3865]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'
[2262]263          CALL abort_gcm('infotrac_init','Bad number of water phases',1)
[3865]264          ENDIF
[3872]265       ENDIF
[1563]266       ! nbtr has been read from INCA by init_const_lmdz() in gcm.F
[2372]267#ifdef INCA
[3870]268       CALL Init_chem_inca_trac(nqINCA)
269#else
270       nqINCA=0
271#endif
272       nbtr=nqINCA+nqCO2
[3872]273       nqtrue=nbtr+nqo
[3865]274       WRITE(lunout,*) trim(modname),': nqo = ',nqo
275       WRITE(lunout,*) trim(modname),': nbtr = ',nbtr
276       WRITE(lunout,*) trim(modname),': nqtrue = ',nqtrue
[3870]277       WRITE(lunout,*) trim(modname),': nqCO2 = ',nqCO2
[3872]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))
[3870]280    ENDIF   ! type_trac 'inca' ou 'inco'
[2262]281!>jyg
[1114]282
[1454]283    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
[3872]284       WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowed. 2 tracers is the minimum'
[1114]285       CALL abort_gcm('infotrac_init','Not enough tracers',1)
[3872]286    ENDIF
[1563]287   
[2262]288!jyg<
[1565]289! Transfert number of tracers to Reprobus
[2262]290!!    IF (type_trac == 'repr') THEN
291!!#ifdef REPROBUS
292!!       CALL Init_chem_rep_trac(nbtr)
293!!#endif
[3872]294!!    ENDIF
[2262]295!>jyg
[1563]296       
[1114]297!
[2262]298! Allocate variables depending on nqtrue
[1114]299!
[2270]300    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
[2566]301
[2262]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
[1114]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)
[3923]321!     iadv = 17   schema  PPM Semi Monotone (overshoots autoris?s)
322!     iadv = 18   schema  PPM Positif Defini (overshoots undershoots autoris?s)
[1114]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!---------------------------------------------------------------------
[3361]335    IF (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac == 'coag' .OR. type_trac == 'co2i') THEN
[1114]336       IF(ierr.EQ.0) THEN
337          ! Continue to read tracer.def
338          DO iq=1,nqtrue
[2270]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.
[2391]352             continu=.true.
353             nouveau_traceurdef=.false.
[2270]354             iiq=1
355             do while (continu)
356                if (tchaine(iiq:iiq).eq.' ') then
[2391]357                  nouveau_traceurdef=.true.
358                  continu=.false.
[2270]359                else if (iiq.lt.LEN_TRIM(tchaine)) then
360                  iiq=iiq+1
361                else
[2391]362                  continu=.false.
[2270]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)
[3923]369                tnom_transp(iq)=tchaine(iiq+1:30)
[2270]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
[3872]379          ENDDO!DO iq=1,nqtrue
[1114]380          CLOSE(90) 
[2270]381
[1454]382       ELSE ! Without tracer.def, set default values
383         if (planet_type=="earth") then
384          ! for Earth, default is to have 4 tracers
[1114]385          hadv(1) = 14
386          vadv(1) = 14
387          tnom_0(1) = 'H2Ov'
[2270]388          tnom_transp(1) = 'air'
[1114]389          hadv(2) = 10
390          vadv(2) = 10
391          tnom_0(2) = 'H2Ol'
[2270]392          tnom_transp(2) = 'air'
[1114]393          hadv(3) = 10
394          vadv(3) = 10
395          tnom_0(3) = 'RN'
[2270]396          tnom_transp(3) = 'air'
[1114]397          hadv(4) = 10
398          vadv(4) = 10
399          tnom_0(4) = 'PB'
[2270]400          tnom_transp(4) = 'air'
[1454]401         else ! default for other planets
402          hadv(1) = 10
403          vadv(1) = 10
404          tnom_0(1) = 'dummy'
[2270]405          tnom_transp(1) = 'dummy'
[1454]406         endif ! of if (planet_type=="earth")
[3872]407       ENDIF
[1114]408       
[1454]409       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
[3865]410       WRITE(lunout,*) trim(modname),': nombre total de traceurs ',nqtrue
[1114]411       DO iq=1,nqtrue
[2270]412          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
[1114]413       END DO
414
[2690]415       IF ( planet_type=='earth') THEN
[2262]416         !CR: nombre de traceurs de l eau
[2690]417         IF (tnom_0(3) == 'H2Oi') THEN
[2262]418            nqo=3
[2690]419         ELSE
[2262]420            nqo=2
[2690]421         ENDIF
[2262]422         ! For Earth, water vapour & liquid tracers are not in the physics
423         nbtr=nqtrue-nqo
[2690]424       ELSE
[2262]425         ! Other planets (for now); we have the same number of tracers
426         ! in the dynamics than in the physics
427         nbtr=nqtrue
[2690]428       ENDIF
[2262]429
[2690]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
[3865]468    ENDIF  ! (type_trac == 'lmdz' .OR. type_trac == 'repr' .OR. type_trac = 'coag' .OR. type_trac = 'co2i')
[2262]469!jyg<
470!
[2566]471
[1114]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
[1454]488          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
[1114]489          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
[3872]490       ENDIF
[1114]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
[1454]498       WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
[1114]499       WRITE(lunout,*) 'makes it necessary to add tracers'
[1454]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'
[1114]502
503    ELSE
504       ! The true number of tracers is also the total number
505       nqtot = nqtrue
[3872]506    ENDIF
[1114]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
[1454]529          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
530
[1114]531          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
[3872]532       ENDIF
[1114]533     
534       str1=tnom_0(iq)
535       tname(new_iq)= tnom_0(iq)
536       IF (iadv(new_iq)==0) THEN
[1454]537          ttext(new_iq)=trim(str1)
[1114]538       ELSE
[1454]539          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
[3872]540       ENDIF
[1114]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
[1454]548             ttext(new_iq)=trim(str2)//txts(jq)
549             tname(new_iq)=trim(str1)//txts(jq)
[1114]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
[1454]555             ttext(new_iq)=trim(str2)//txtp(jq)
556             tname(new_iq)=trim(str1)//txtp(jq)
[1114]557          END DO
[3872]558       ENDIF
[1114]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
[1454]575    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
576    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
[1114]577    DO iq=1,nqtot
[1454]578       WRITE(lunout,*) iadv(iq),niadv(iq),&
579       ' ',trim(tname(iq)),' ',trim(ttext(iq))
[1114]580    END DO
581
[1279]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
[1454]588          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
[1279]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
[1454]591          WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
[1279]592          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
[3872]593       ENDIF
[1279]594    END DO
595
[2270]596
[3923]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
[2270]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
[3923]611        ! ceci est un traceur p?re
[2270]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
[3923]616        ! ceci est un fils. Qui est son p?re?
[2270]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
[3923]622            ! Son p?re est ipere
[2270]623            WRITE(lunout,*) 'Le traceur',iq,'appele ', &
624      &          trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
[3923]625            if (iq.eq.ipere) then
626                CALL abort_gcm('infotrac_init','Un fils est son propre pere',1)
627            endif
[2270]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 ', &
[3361]636      &          trim(tnom_0(iq)),', est orphelin.'
[2270]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
[3923]648! Calculer le nombre de descendants ? partir de iqfils et de nbfils
[2270]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
[3923]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
[2270]673    do iq=1,nqtot
674      if (iqpere(iq).gt.0) then
[3923]675        ! ce traceur a un p?re qui n'est pas l'air
676        ! Seul le sch?ma 10 est autoris?
[2270]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
[3923]681        ! Le traceur p?re ne peut ?tre advect? que par sch?ma 10 ou 14:
[2270]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
[3923]690
[2270]691! detecter quels sont les traceurs isotopiques parmi des traceurs
692    call infotrac_isoinit(tnom_0,nqtrue)
[3923]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
[1114]827!-----------------------------------------------------------------------
828! Finalize :
829!
[2270]830    DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
[1114]831
[3923]832    WRITE(lunout,*) 'infotrac init fin'
[2180]833
[1114]834  END SUBROUTINE infotrac_init
835
[2270]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
[3923]848    CHARACTER(len=30) tnom_0(nqtrue)
[2270]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
[2372]883!       write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
[2270]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))
[2372]888!         write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac     
[2270]889         IF (tnom_0(iq) == tnom_trac) then
[2372]890!          write(lunout,*) 'Ce traceur est un isotope'
[2270]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
[2372]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)
[2270]901          goto 20
902         else if (iqpere(iq).gt.0) then         
903          if (tnom_0(iqpere(iq)) == tnom_trac) then
[2372]904!           write(lunout,*) 'Ce traceur est le fils d''un isotope'
[2270]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
[2372]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)
[2270]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
[2372]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
[2270]927
928      do ixt= 1,niso_possibles 
929
930        if (nb_iso(ixt,1).eq.1) then
[3923]931          ! on v?rifie que toutes les phases ont le m?me nombre de
[2270]932          ! traceurs
933          do phase=2,nqo
934            if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
[2372]935!              write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
[2270]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
[3923]944          ! on v?rifie que toutes les phases ont le m?me nombre de
[2270]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
[3923]953          ! on v?rifie que tous les isotopes ont le m?me nombre de
[2270]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)
[2372]974!    WRITE(lunout,*) 'niso=',niso
975!    WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso   
[2270]976 
977    ! flags isotopiques:
978    if (niso.gt.0) then
979        ok_isotopes=.true.
980    else
981        ok_isotopes=.false.
982    endif
[2372]983!    WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
[2270]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 
[2372]993!    WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
994!    WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
[2270]995
996    if (ntraceurs_zone.gt.0) then
997        ok_isotrac=.true.
998    else
999        ok_isotrac=.false.
1000    endif   
[2372]1001!    WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
[2270]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
[2372]1012!    WRITE(lunout,*) 'iqiso=',iqiso
[2270]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
[2372]1025!    write(lunout,*) 'index_trac=',index_trac   
[2270]1026
1027! Finalize :
1028    DEALLOCATE(nb_iso)
1029
1030  END SUBROUTINE infotrac_isoinit
1031
[1114]1032END MODULE infotrac
Note: See TracBrowser for help on using the repository browser.