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

Last change on this file since 4046 was 4046, checked in by dcugnet, 2 years ago

First commit for new tracers.

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