source: trunk/LMDZ.COMMON/libf/dyn3d_common/infotrac.F90 @ 3000

Last change on this file since 3000 was 2604, checked in by emillour, 3 years ago

Common dynamics:
Update "infotrac_init": stop with an error message if file traceur.def is
not found (rather than make the assumption that the model is running
with a dummy tracer)
EM

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