source: LMDZ5/trunk/libf/dyn3d_common/infotrac.F90 @ 3356

Last change on this file since 3356 was 2690, checked in by oboucher, 8 years ago

Adding a module for stratospheric aerosols with a bin scheme.
The module gets activated with -strataer true compiling option.
May not quite work yet, more testing needed, but should not affect
the rest of LMDz as everything is under a CPP_StratAer key.

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