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

Last change on this file since 2567 was 2567, checked in by acozic, 8 years ago

correct a bug / hadv_inca and vadv_inca need to be defined even when we don't use inca model

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