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

Last change on this file since 3639 was 3362, checked in by oboucher, 6 years ago

removing two write statements

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