Ignore:
Timestamp:
Jun 14, 2015, 9:13:32 PM (9 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes -r2237:2291 into testing branch

Location:
LMDZ5/branches/testing
Files:
11 deleted
3 edited
7 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dyn3d_common/diagedyn.F

    r1999 r2298  
    5353c======================================================================
    5454 
     55      USE control_mod, ONLY : planet_type
     56     
    5557      IMPLICIT NONE
    5658C
    … …  
    6062#include "iniprint.h"
    6163
    62 #ifdef CPP_EARTH
    63 #include "../phylmd/YOMCST.h"
    64 #include "../phylmd/YOETHF.h"
    65 #endif
     64!#ifdef CPP_EARTH
     65!#include "../phylmd/YOMCST.h"
     66!#include "../phylmd/YOETHF.h"
     67!#endif
     68! Ehouarn: for now set these parameters to what is in Earth physics...
     69!          (cf ../phylmd/suphel.h)
     70!          this should be generalized...
     71      REAL,PARAMETER :: RCPD=
     72     &               3.5*(1000.*(6.0221367E+23*1.380658E-23)/28.9644)
     73      REAL,PARAMETER :: RCPV=
     74     &               4.*(1000.*(6.0221367E+23*1.380658E-23)/18.0153)
     75      REAL,PARAMETER :: RCS=RCPV
     76      REAL,PARAMETER :: RCW=RCPV
     77      REAL,PARAMETER :: RLSTT=2.8345E+6
     78      REAL,PARAMETER :: RLVTT=2.5008E+6
     79!
    6680C
    6781      INTEGER imjmp1
    … …  
    140154
    141155
    142 #ifdef CPP_EARTH
     156!#ifdef CPP_EARTH
     157      IF (planet_type=="earth") THEN
     158     
    143159c======================================================================
    144160C     Compute Kinetic enrgy
    … …  
    314330      ec_pre (idiag)    = ec_tot
    315331C
    316 #else
    317       write(lunout,*)'diagedyn: Needs Earth physics to function'
    318 #endif
     332!#else
     333      ELSE
     334        write(lunout,*)'diagedyn: set to function with Earth parameters'
     335      ENDIF ! of if (planet_type=="earth")
     336!#endif
    319337! #endif of #ifdef CPP_EARTH
    320338      RETURN

  • LMDZ5/branches/testing/libf/dyn3d_common/infotrac.F90

    r2187 r2298  
    1212  INTEGER, SAVE :: nbtr
    1313
     14! CRisi: nb traceurs pères= directement advectés par l'air
     15  INTEGER, SAVE :: nqperes
     16
    1417! Name variables
    1518  CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics
    … …  
    2225!         dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code.
    2326  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
    2434
    2535! conv_flg(it)=0 : convection desactivated for tracer number it
    … …  
    3040  CHARACTER(len=4),SAVE :: type_trac
    3141  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
    3257 
    3358CONTAINS
    … …  
    6388
    6489    CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_0  ! tracer short name
     90    CHARACTER(len=15), ALLOCATABLE, DIMENSION(:) :: tnom_transp ! transporting fluid short name: CRisi
    6591    CHARACTER(len=3), DIMENSION(30) :: descrq
    6692    CHARACTER(len=1), DIMENSION(3)  :: txts
    … …  
    7096    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
    7197    INTEGER :: iq, new_iq, iiq, jq, ierr
     98    INTEGER :: ifils,ipere,generation ! CRisi
     99    LOGICAL :: continu,nouveau_traceurdef
     100    INTEGER :: IOstatus ! gestion de la retrocompatibilite de traceur.def
     101    CHARACTER(len=15) :: tchaine   
    72102
    73103    character(len=*),parameter :: modname="infotrac_init"
    … …  
    134164          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
    135165          READ(90,*) nqtrue
     166          write(lunout,*) 'nqtrue=',nqtrue
    136167       ELSE
    137168          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
    … …  
    143174          endif
    144175       END IF
    145        if ( planet_type=='earth') then
    146          ! For Earth, water vapour & liquid tracers are not in the physics
    147          nbtr=nqtrue-2
    148        else
    149          ! Other planets (for now); we have the same number of tracers
    150          ! in the dynamics than in the physics
    151          nbtr=nqtrue
    152        endif
     176!jyg<
     177!!       if ( planet_type=='earth') then
     178!!         ! For Earth, water vapour & liquid tracers are not in the physics
     179!!         nbtr=nqtrue-2
     180!!       else
     181!!         ! Other planets (for now); we have the same number of tracers
     182!!         ! in the dynamics than in the physics
     183!!         nbtr=nqtrue
     184!!       endif
     185!>jyg
    153186    ELSE ! type_trac=inca
     187!jyg<
     188       ! The traceur.def file is used to define the number "nqo" of water phases
     189       ! present in the simulation. Default : nqo = 2.
     190       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
     191       IF(ierr.EQ.0) THEN
     192          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
     193          READ(90,*) nqo
     194       ELSE
     195          WRITE(lunout,*) trim(modname),': Using default value for nqo'
     196          nqo=2
     197       ENDIF
     198       IF (nqo /= 2 .OR. nqo /= 3 ) THEN
     199          WRITE(lunout,*) trim(modname),': nqo=',nqo, ' is not allowded. Only 2 or 3 water phases allowed'
     200          CALL abort_gcm('infotrac_init','Bad number of water phases',1)
     201       END IF
    154202       ! nbtr has been read from INCA by init_const_lmdz() in gcm.F
    155        nqtrue=nbtr+2
    156     END IF
     203       nqtrue=nbtr+nqo
     204!!       nqtrue=nbtr+2
     205    END IF   ! type_trac
     206!>jyg
    157207
    158208    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
    … …  
    161211    END IF
    162212   
     213!jyg<
    163214! Transfert number of tracers to Reprobus
    164     IF (type_trac == 'repr') THEN
    165 #ifdef REPROBUS
    166        CALL Init_chem_rep_trac(nbtr)
    167 #endif
    168     END IF
     215!!    IF (type_trac == 'repr') THEN
     216!!#ifdef REPROBUS
     217!!       CALL Init_chem_rep_trac(nbtr)
     218!!#endif
     219!!    END IF
     220!>jyg
    169221       
    170222!
    171 ! Allocate variables depending on nqtrue and nbtr
    172 !
    173     ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue))
    174     ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
    175     conv_flg(:) = 1 ! convection activated for all tracers
    176     pbl_flg(:)  = 1 ! boundary layer activated for all tracers
     223! Allocate variables depending on nqtrue
     224!
     225    ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue),tnom_transp(nqtrue))
     226!
     227!jyg<
     228!!    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
     229!!    conv_flg(:) = 1 ! convection activated for all tracers
     230!!    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
     231!>jyg
    177232
    178233!-----------------------------------------------------------------------
    … …  
    206261          ! Continue to read tracer.def
    207262          DO iq=1,nqtrue
    208              READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    209           END DO
     263
     264             write(*,*) 'infotrac 237: iq=',iq
     265             ! CRisi: ajout du nom du fluide transporteur
     266             ! mais rester retro compatible
     267             READ(90,'(I2,X,I2,X,A)',IOSTAT=IOstatus) hadv(iq),vadv(iq),tchaine
     268             write(lunout,*) 'iq,hadv(iq),vadv(iq)=',iq,hadv(iq),vadv(iq)
     269             write(lunout,*) 'tchaine=',trim(tchaine)
     270             write(*,*) 'infotrac 238: IOstatus=',IOstatus
     271             if (IOstatus.ne.0) then
     272                CALL abort_gcm('infotrac_init','Pb dans la lecture de traceur.def',1)
     273             endif
     274             ! Y-a-t-il 1 ou 2 noms de traceurs? -> On regarde s'il y a un
     275             ! espace ou pas au milieu de la chaine.
     276             continu=1
     277             nouveau_traceurdef=0
     278             iiq=1
     279             do while (continu)
     280                if (tchaine(iiq:iiq).eq.' ') then
     281                  nouveau_traceurdef=1
     282                  continu=0
     283                else if (iiq.lt.LEN_TRIM(tchaine)) then
     284                  iiq=iiq+1
     285                else
     286                  continu=0     
     287                endif
     288             enddo
     289             write(*,*) 'iiq,nouveau_traceurdef=',iiq,nouveau_traceurdef
     290             if (nouveau_traceurdef) then
     291                write(lunout,*) 'C''est la nouvelle version de traceur.def'
     292                tnom_0(iq)=tchaine(1:iiq-1)
     293                tnom_transp(iq)=tchaine(iiq+1:15)
     294             else
     295                write(lunout,*) 'C''est l''ancienne version de traceur.def'
     296                write(lunout,*) 'On suppose que les traceurs sont tous d''air'
     297                tnom_0(iq)=tchaine
     298                tnom_transp(iq) = 'air'
     299             endif
     300             write(lunout,*) 'tnom_0(iq)=<',trim(tnom_0(iq)),'>'
     301             write(lunout,*) 'tnom_transp(iq)=<',trim(tnom_transp(iq)),'>'
     302
     303          END DO !DO iq=1,nqtrue
    210304          CLOSE(90) 
     305
    211306       ELSE ! Without tracer.def, set default values
    212307         if (planet_type=="earth") then
    … …  
    215310          vadv(1) = 14
    216311          tnom_0(1) = 'H2Ov'
     312          tnom_transp(1) = 'air'
    217313          hadv(2) = 10
    218314          vadv(2) = 10
    219315          tnom_0(2) = 'H2Ol'
     316          tnom_transp(2) = 'air'
    220317          hadv(3) = 10
    221318          vadv(3) = 10
    222319          tnom_0(3) = 'RN'
     320          tnom_transp(3) = 'air'
    223321          hadv(4) = 10
    224322          vadv(4) = 10
    225323          tnom_0(4) = 'PB'
     324          tnom_transp(4) = 'air'
    226325         else ! default for other planets
    227326          hadv(1) = 10
    228327          vadv(1) = 10
    229328          tnom_0(1) = 'dummy'
     329          tnom_transp(1) = 'dummy'
    230330         endif ! of if (planet_type=="earth")
    231331       END IF
    232 
    233 !CR: nombre de traceurs de l eau
    234        if (tnom_0(3) == 'H2Oi') then
    235           nqo=3
    236        else
    237           nqo=2
    238        endif
    239332       
    240333       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
    241334       WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
    242335       DO iq=1,nqtrue
    243           WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
     336          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
    244337       END DO
    245338
    246     ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
     339       if ( planet_type=='earth') then
     340         !CR: nombre de traceurs de l eau
     341         if (tnom_0(3) == 'H2Oi') then
     342            nqo=3
     343         else
     344            nqo=2
     345         endif
     346         ! For Earth, water vapour & liquid tracers are not in the physics
     347         nbtr=nqtrue-nqo
     348       else
     349         ! Other planets (for now); we have the same number of tracers
     350         ! in the dynamics than in the physics
     351         nbtr=nqtrue
     352       endif
     353
     354    ENDIF  ! (type_trac == 'lmdz' .OR. type_trac == 'repr')
     355!jyg<
     356!
     357! Transfert number of tracers to Reprobus
     358    IF (type_trac == 'repr') THEN
     359#ifdef REPROBUS
     360       CALL Init_chem_rep_trac(nbtr)
     361#endif
     362    END IF
     363!
     364! Allocate variables depending on nbtr
     365!
     366    ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
     367    conv_flg(:) = 1 ! convection activated for all tracers
     368    pbl_flg(:)  = 1 ! boundary layer activated for all tracers
     369!
     370!!    ELSE  ! type_trac=inca : config_inca='aero' ou 'chem'
     371!
     372    IF (type_trac == 'inca') THEN   ! config_inca='aero' ou 'chem'
     373!>jyg
    247374! le module de chimie fournit les noms des traceurs
    248375! et les schemas d'advection associes.
    … …  
    258385       tnom_0(1)='H2Ov'
    259386       tnom_0(2)='H2Ol'
    260 
    261        DO iq =3,nqtrue
    262           tnom_0(iq)=solsym(iq-2)
     387       IF (nqo == 3) tnom_0(3)='H2Oi'     !! jyg
     388
     389!jyg<
     390       DO iq = nqo+1, nqtrue
     391          tnom_0(iq)=solsym(iq-nqo)
    263392       END DO
    264        nqo = 2
    265 
    266     END IF ! type_trac
     393!!       DO iq =3,nqtrue
     394!!          tnom_0(iq)=solsym(iq-2)
     395!!       END DO
     396!!       nqo = 2
     397!>jyg
     398
     399    END IF ! (type_trac == 'inca')
    267400
    268401!-----------------------------------------------------------------------
    … …  
    390523    END DO
    391524
     525
     526! CRisi: quels sont les traceurs fils et les traceurs pères.
     527! initialiser tous les tableaux d'indices liés aux traceurs familiaux
     528! + vérifier que tous les pères sont écrits en premières positions
     529    ALLOCATE(nqfils(nqtot),nqdesc(nqtot))   
     530    ALLOCATE(iqfils(nqtot,nqtot))   
     531    ALLOCATE(iqpere(nqtot))
     532    nqperes=0
     533    nqfils(:)=0
     534    nqdesc(:)=0
     535    iqfils(:,:)=0
     536    iqpere(:)=0
     537    nqdesc_tot=0   
     538    DO iq=1,nqtot
     539      if (tnom_transp(iq) == 'air') then
     540        ! ceci est un traceur père
     541        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un pere'
     542        nqperes=nqperes+1
     543        iqpere(iq)=0
     544      else !if (tnom_transp(iq) == 'air') then
     545        ! ceci est un fils. Qui est son père?
     546        WRITE(lunout,*) 'Le traceur',iq,', appele ',trim(tnom_0(iq)),', est un fils'
     547        continu=.true.
     548        ipere=1
     549        do while (continu)           
     550          if (tnom_transp(iq) == tnom_0(ipere)) then
     551            ! Son père est ipere
     552            WRITE(lunout,*) 'Le traceur',iq,'appele ', &
     553      &          trim(tnom_0(iq)),' est le fils de ',ipere,'appele ',trim(tnom_0(ipere))
     554            nqfils(ipere)=nqfils(ipere)+1 
     555            iqfils(nqfils(ipere),ipere)=iq
     556            iqpere(iq)=ipere         
     557            continu=.false.
     558          else !if (tnom_transp(iq) == tnom_0(ipere)) then
     559            ipere=ipere+1
     560            if (ipere.gt.nqtot) then
     561                WRITE(lunout,*) 'Le traceur',iq,'appele ', &
     562      &          trim(tnom_0(iq)),', est orpelin.'
     563                CALL abort_gcm('infotrac_init','Un traceur est orphelin',1)
     564            endif !if (ipere.gt.nqtot) then
     565          endif !if (tnom_transp(iq) == tnom_0(ipere)) then
     566        enddo !do while (continu)
     567      endif !if (tnom_transp(iq) == 'air') then
     568    enddo !DO iq=1,nqtot
     569    WRITE(lunout,*) 'infotrac: nqperes=',nqperes   
     570    WRITE(lunout,*) 'nqfils=',nqfils
     571    WRITE(lunout,*) 'iqpere=',iqpere
     572    WRITE(lunout,*) 'iqfils=',iqfils
     573
     574! Calculer le nombre de descendants à partir de iqfils et de nbfils
     575    DO iq=1,nqtot   
     576      generation=0
     577      continu=.true.
     578      ifils=iq
     579      do while (continu)
     580        ipere=iqpere(ifils)
     581        if (ipere.gt.0) then
     582         nqdesc(ipere)=nqdesc(ipere)+1   
     583         nqdesc_tot=nqdesc_tot+1     
     584         iqfils(nqdesc(ipere),ipere)=iq
     585         ifils=ipere
     586         generation=generation+1
     587        else !if (ipere.gt.0) then
     588         continu=.false.
     589        endif !if (ipere.gt.0) then
     590      enddo !do while (continu)   
     591      WRITE(lunout,*) 'Le traceur ',iq,', appele ',trim(tnom_0(iq)),' est un traceur de generation: ',generation
     592    enddo !DO iq=1,nqtot
     593    WRITE(lunout,*) 'infotrac: nqdesc=',nqdesc
     594    WRITE(lunout,*) 'iqfils=',iqfils
     595    WRITE(lunout,*) 'nqdesc_tot=',nqdesc_tot
     596
     597! Interdire autres schémas que 10 pour les traceurs fils, et autres schémas
     598! que 10 et 14 si des pères ont des fils
     599    do iq=1,nqtot
     600      if (iqpere(iq).gt.0) then
     601        ! ce traceur a un père qui n'est pas l'air
     602        ! Seul le schéma 10 est autorisé
     603        if (iadv(iq)/=10) then
     604           WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for sons'
     605          CALL abort_gcm('infotrac_init','Sons should be advected by scheme 10',1)
     606        endif
     607        ! Le traceur père ne peut être advecté que par schéma 10 ou 14:
     608        IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
     609          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not implemented for fathers'
     610          CALL abort_gcm('infotrac_init','Fathers should be advected by scheme 10 ou 14',1)
     611        endif !IF (iadv(iqpere(iq))/=10 .AND. iadv(iqpere(iq))/=14) THEN
     612     endif !if (iqpere(iq).gt.0) the
     613    enddo !do iq=1,nqtot
     614
     615
     616! detecter quels sont les traceurs isotopiques parmi des traceurs
     617    call infotrac_isoinit(tnom_0,nqtrue)
     618       
    392619!-----------------------------------------------------------------------
    393620! Finalize :
    394621!
    395     DEALLOCATE(tnom_0, hadv, vadv)
     622    DEALLOCATE(tnom_0, hadv, vadv,tnom_transp)
    396623
    397624
    398625  END SUBROUTINE infotrac_init
    399626
     627  SUBROUTINE infotrac_isoinit(tnom_0,nqtrue)
     628
     629#ifdef CPP_IOIPSL
     630  use IOIPSL
     631#else
     632  ! if not using IOIPSL, we still need to use (a local version of) getin
     633  use ioipsl_getincom
     634#endif
     635  implicit none
     636 
     637    ! inputs
     638    INTEGER nqtrue
     639    CHARACTER(len=15) tnom_0(nqtrue)
     640   
     641    ! locals   
     642    CHARACTER(len=3), DIMENSION(niso_possibles) :: tnom_iso
     643    INTEGER, ALLOCATABLE,DIMENSION(:,:) :: nb_iso,nb_traciso
     644    INTEGER, ALLOCATABLE,DIMENSION(:) :: nb_isoind
     645    INTEGER :: ntraceurs_zone_prec,iq,phase,ixt,iiso,izone
     646    CHARACTER(len=19) :: tnom_trac
     647    INCLUDE "iniprint.h"
     648
     649    tnom_iso=(/'eau','HDO','O18','O17','HTO'/)
     650
     651    ALLOCATE(nb_iso(niso_possibles,nqo))
     652    ALLOCATE(nb_isoind(nqo))
     653    ALLOCATE(nb_traciso(niso_possibles,nqo))
     654    ALLOCATE(iso_num(nqtot))
     655    ALLOCATE(iso_indnum(nqtot))
     656    ALLOCATE(zone_num(nqtot))
     657    ALLOCATE(phase_num(nqtot))
     658     
     659    iso_num(:)=0
     660    iso_indnum(:)=0
     661    zone_num(:)=0
     662    phase_num(:)=0
     663    indnum_fn_num(:)=0
     664    use_iso(:)=.false. 
     665    nb_iso(:,:)=0 
     666    nb_isoind(:)=0     
     667    nb_traciso(:,:)=0
     668    niso=0
     669    ntraceurs_zone=0 
     670    ntraceurs_zone_prec=0
     671    ntraciso=0
     672
     673    do iq=nqo+1,nqtot
     674       write(lunout,*) 'infotrac 569: iq,tnom_0(iq)=',iq,tnom_0(iq)
     675       do phase=1,nqo   
     676        do ixt= 1,niso_possibles   
     677         tnom_trac=trim(tnom_0(phase))//'_'
     678         tnom_trac=trim(tnom_trac)//trim(tnom_iso(ixt))
     679         write(*,*) 'phase,ixt,tnom_trac=',phase,ixt,tnom_trac     
     680         IF (tnom_0(iq) == tnom_trac) then
     681          write(lunout,*) 'Ce traceur est un isotope'
     682          nb_iso(ixt,phase)=nb_iso(ixt,phase)+1   
     683          nb_isoind(phase)=nb_isoind(phase)+1   
     684          iso_num(iq)=ixt
     685          iso_indnum(iq)=nb_isoind(phase)
     686          indnum_fn_num(ixt)=iso_indnum(iq)
     687          phase_num(iq)=phase
     688          write(lunout,*) 'iso_num(iq)=',iso_num(iq)
     689          write(lunout,*) 'iso_indnum(iq)=',iso_indnum(iq)
     690          write(lunout,*) 'indnum_fn_num(ixt)=',indnum_fn_num(ixt)
     691          write(lunout,*) 'phase_num(iq)=',phase_num(iq)
     692          goto 20
     693         else if (iqpere(iq).gt.0) then         
     694          if (tnom_0(iqpere(iq)) == tnom_trac) then
     695           write(lunout,*) 'Ce traceur est le fils d''un isotope'
     696           ! c'est un traceur d'isotope
     697           nb_traciso(ixt,phase)=nb_traciso(ixt,phase)+1
     698           iso_num(iq)=ixt
     699           iso_indnum(iq)=indnum_fn_num(ixt)
     700           zone_num(iq)=nb_traciso(ixt,phase)
     701           phase_num(iq)=phase
     702           write(lunout,*) 'iso_num(iq)=',iso_num(iq)
     703           write(lunout,*) 'phase_num(iq)=',phase_num(iq)
     704           write(lunout,*) 'zone_num(iq)=',zone_num(iq)
     705           goto 20
     706          endif !if (tnom_0(iqpere(iq)) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
     707         endif !IF (tnom_0(iq) == trim(tnom_0(phase))//trim(tnom_iso(ixt))) then
     708        enddo !do ixt= niso_possibles
     709       enddo !do phase=1,nqo
     710  20   continue
     711      enddo !do iq=1,nqtot
     712
     713      write(lunout,*) 'iso_num=',iso_num
     714      write(lunout,*) 'iso_indnum=',iso_indnum
     715      write(lunout,*) 'zone_num=',zone_num 
     716      write(lunout,*) 'phase_num=',phase_num
     717      write(lunout,*) 'indnum_fn_num=',indnum_fn_num
     718
     719      do ixt= 1,niso_possibles 
     720
     721        if (nb_iso(ixt,1).eq.1) then
     722          ! on vérifie que toutes les phases ont le même nombre de
     723          ! traceurs
     724          do phase=2,nqo
     725            if (nb_iso(ixt,phase).ne.nb_iso(ixt,1)) then
     726              write(lunout,*) 'ixt,phase,nb_iso=',ixt,phase,nb_iso(ixt,phase)
     727              CALL abort_gcm('infotrac_init','Phases must have same number of isotopes',1)
     728            endif
     729          enddo !do phase=2,nqo
     730
     731          niso=niso+1
     732          use_iso(ixt)=.true.
     733          ntraceurs_zone=nb_traciso(ixt,1)
     734
     735          ! on vérifie que toutes les phases ont le même nombre de
     736          ! traceurs
     737          do phase=2,nqo
     738            if (nb_traciso(ixt,phase).ne.ntraceurs_zone) then
     739              write(lunout,*) 'ixt,phase,nb_traciso=',ixt,phase,nb_traciso(ixt,phase)
     740              write(lunout,*) 'ntraceurs_zone=',ntraceurs_zone
     741              CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
     742            endif 
     743          enddo  !do phase=2,nqo
     744          ! on vérifie que tous les isotopes ont le même nombre de
     745          ! traceurs
     746          if (ntraceurs_zone_prec.gt.0) then               
     747            if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
     748              ntraceurs_zone_prec=ntraceurs_zone
     749            else !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
     750              write(*,*) 'ntraceurs_zone_prec,ntraceurs_zone=',ntraceurs_zone_prec,ntraceurs_zone   
     751              CALL abort_gcm('infotrac_init', &
     752               &'Isotope tracers are not well defined in traceur.def',1)           
     753            endif !if (ntraceurs_zone.eq.ntraceurs_zone_prec) then
     754           endif !if (ntraceurs_zone_prec.gt.0) then
     755
     756        else if (nb_iso(ixt,1).ne.0) then
     757           WRITE(lunout,*) 'nqo,ixt=',nqo,ixt
     758           WRITE(lunout,*) 'nb_iso(ixt,1)=',nb_iso(ixt,1)   
     759           CALL abort_gcm('infotrac_init','Isotopes are not well defined in traceur.def',1)     
     760        endif   !if (nb_iso(ixt,1).eq.1) then       
     761    enddo ! do ixt= niso_possibles
     762
     763    ! dimensions isotopique:
     764    ntraciso=niso*(ntraceurs_zone+1)
     765    WRITE(lunout,*) 'niso=',niso
     766    WRITE(lunout,*) 'ntraceurs_zone,ntraciso=',ntraceurs_zone,ntraciso   
     767 
     768    ! flags isotopiques:
     769    if (niso.gt.0) then
     770        ok_isotopes=.true.
     771    else
     772        ok_isotopes=.false.
     773    endif
     774    WRITE(lunout,*) 'ok_isotopes=',ok_isotopes
     775 
     776    if (ok_isotopes) then
     777        ok_iso_verif=.false.
     778        call getin('ok_iso_verif',ok_iso_verif)
     779        ok_init_iso=.false.
     780        call getin('ok_init_iso',ok_init_iso)
     781        tnat=(/1.0,155.76e-6,2005.2e-6,0.004/100.,0.0/)
     782        alpha_ideal=(/1.0,1.01,1.006,1.003,1.0/)
     783    endif !if (ok_isotopes) then 
     784    WRITE(lunout,*) 'ok_iso_verif=',ok_iso_verif
     785    WRITE(lunout,*) 'ok_init_iso=',ok_init_iso
     786
     787    if (ntraceurs_zone.gt.0) then
     788        ok_isotrac=.true.
     789    else
     790        ok_isotrac=.false.
     791    endif   
     792    WRITE(lunout,*) 'ok_isotrac=',ok_isotrac
     793
     794    ! remplissage du tableau iqiso(ntraciso,phase)
     795    ALLOCATE(iqiso(ntraciso,nqo))   
     796    iqiso(:,:)=0     
     797    do iq=1,nqtot
     798        if (iso_num(iq).gt.0) then
     799          ixt=iso_indnum(iq)+zone_num(iq)*niso
     800          iqiso(ixt,phase_num(iq))=iq
     801        endif
     802    enddo
     803    WRITE(lunout,*) 'iqiso=',iqiso
     804
     805    ! replissage du tableau index_trac(ntraceurs_zone,niso)
     806    ALLOCATE(index_trac(ntraceurs_zone,niso)) 
     807    if (ok_isotrac) then
     808        do iiso=1,niso
     809          do izone=1,ntraceurs_zone
     810             index_trac(izone,iiso)=iiso+izone*niso
     811          enddo
     812        enddo
     813    else !if (ok_isotrac) then     
     814        index_trac(:,:)=0.0
     815    endif !if (ok_isotrac) then
     816    write(lunout,*) 'index_trac=',index_trac   
     817
     818! Finalize :
     819    DEALLOCATE(nb_iso)
     820
     821  END SUBROUTINE infotrac_isoinit
     822
    400823END MODULE infotrac
Note: See TracChangeset for help on using the changeset viewer.