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

Last change on this file since 2048 was 1971, checked in by emillour, 7 years ago

Common dynamics:
Enable tracer names to be up to 30 characters long (instead of 20).
EM

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