Ignore:
Timestamp:
Jan 29, 2018, 11:41:36 PM (7 years ago)
Author:
jvatant
Message:

+ Generalize the use of nkim parameter in comchem_h instead
of having x=44 hardcoded everywhere in the code.
If ever one day we have nkim!=44 it will be useful !
+ Also add the other parameters for the chem. in comchem_h
+ Some rearrangment in physiq_mod of chemistry variables

-> especially less unsane "nq-nmicro" -> nkim instead.

+ In calchim, preliminary modifs to adapt to the new managment
of i/o chemistry implemented in last commits.
--JVO

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/calchim.F90

    r1795 r1903  
    1       SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad, &
    2                          dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc,NLEV)
     1      SUBROUTINE calchim(ngrid,qy_c,declin_rad,ls_rad, &
     2                         dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc)
    33     
    44!--------------------------------------------------------------------
    … …  
    1717!
    1818     
     19      use comchem_h
    1920      use dimphy
    20      
    2121      use datafile_mod, only: datadir
    22      
    2322      use comcstfi_mod, only: rad, pi, kbol
    2423      use moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
    … …  
    2726      implicit none
    2827
    29 !   Parameters ( dans common_mod in old titan)
    30 !              + doivent etre en accord avec titan.h
    31 !   -------------------------------------------------
    32 
    33       INTEGER,PARAMETER ::   NC   = 44     ! Nb de composes dans la chimie
    34       INTEGER,PARAMETER ::   ND   = 54     ! Nb de photodissociations
    35       INTEGER,PARAMETER ::   NR   = 377    ! Nb de reactions dans la chimie
    36       INTEGER,PARAMETER ::   NLRT = 650    ! Pour l'UV 650 niveaux de 2 km
    37      
    38 !     Dummy parameters without microphysics
    39 !      + Just here to keep the whole stuff running without modif C sources
    40 !     ---------------------------------------------------------------------
    41 
    42       INTEGER  :: utilaer(16)
    43       INTEGER  :: aerprod = 0
    44       INTEGER  :: htoh2   = 0
    45 
     28! ------------------------------------------
     29! *********** 0. Declarations *************
     30! ------------------------------------------
     31     
    4632!    Arguments
    4733!    ---------
    4834
    49       INTEGER      nlon                   ! nb of horiz points
    50       INTEGER      ny                     ! nb de composes (nqmax-nmicro)
    51       REAL*8         qy_c(nlon,klev,NC)     ! Especes chimiques apres adv.+diss.
    52       character*10 nomqy_c(NC+1)          ! Noms des especes chimiques
     35      INTEGER      ngrid                   ! nb of horiz points
     36      REAL*8         qy_c(ngrid,klev,nkim)     ! Especes chimiques apres adv.+diss.
    5337      REAL*8         declin_rad,ls_rad      ! declinaison et long solaire en radians
    5438      REAL*8         dtchim                 ! pas de temps chimie
    55       REAL*8         ctemp(nlon,klev)       ! Temperature
    56       REAL*8         cplay(nlon,klev)       ! pression (Pa)
    57       REAL*8         cplev(nlon,klev+1)     ! pression intercouches (Pa)
    58       REAL*8         czlay(nlon,klev)       ! altitude (m)
    59       REAL*8         czlev(nlon,klev+1)     ! altitude intercouches (m)
    60      
    61       REAL*8         dqyc(nlon,klev,NC)     ! Tendances especes chimiques
    62      
    63       INTEGER      NLEV                   ! nbp_lev+70 - a changer si =/ 55 layers ?
     39      REAL*8         ctemp(ngrid,klev)       ! Temperature
     40      REAL*8         cplay(ngrid,klev)       ! pression (Pa)
     41      REAL*8         cplev(ngrid,klev+1)     ! pression intercouches (Pa)
     42      REAL*8         czlay(ngrid,klev)       ! altitude (m)
     43      REAL*8         czlev(ngrid,klev+1)     ! altitude intercouches (m)
     44     
     45      REAL*8         dqyc(ngrid,klev,nkim)     ! Tendances especes chimiques
     46     
    6447     
    6548!    Local variables :
    … …  
    7053! variables envoyees dans la chimie: double precision
    7154
    72       REAL*8  temp_c(NLEV)
    73       REAL*8  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
    74       REAL*8  cqy(NLEV,NC)    ! frac mol qui seront modifiees
    75       REAL*8  cqy0(NLEV,NC)    ! frac mol avant modif
    76       REAL*8  surfhaze(NLEV)
    77       REAL*8  cprodaer(NLEV,4),cmaer(NLEV,4)
    78       REAL*8  ccsn(NLEV,4),ccsh(NLEV,4)
     55      REAL*8  temp_c(nlaykim_tot)
     56      REAL*8  press_c(nlaykim_tot)   ! T,p(mbar) a 1 lat donnee
     57      REAL*8  cqy(nlaykim_tot,nkim)    ! frac mol qui seront modifiees
     58      REAL*8  cqy0(nlaykim_tot,nkim)    ! frac mol avant modif
     59      REAL*8  surfhaze(nlaykim_tot)
     60      REAL*8  cprodaer(nlaykim_tot,4),cmaer(nlaykim_tot,4)
     61      REAL*8  ccsn(nlaykim_tot,4),ccsh(nlaykim_tot,4)
    7962! rmil: milieu de couche, grille pour K,D,p,ct,T,y
    8063! rinter: intercouche (grille RA dans la chimie)
    81       REAL*8  rmil(NLEV),rinter(NLEV),nb(NLEV)
     64      REAL*8  rmil(nlaykim_tot),rinter(nlaykim_tot),nb(nlaykim_tot)
    8265      REAL*8,save,allocatable :: kedd(:)
    8366
    84       REAL*8  a,b
    8567      character str1*1,str2*2
    86       REAL*8  latit
    87       character*20 formt1,formt2,formt3
    8868     
    8969!     variables locales initialisees au premier appel
    … …  
    9777      real*8  factalt,factdec,krpddec,krpddecp1,krpddecm1
    9878      real*8  duree
    99       REAL*8,save :: mass(NC)
     79      REAL*8,save :: mass(nkim)
    10080      REAL*8,save,allocatable :: md(:,:)
    10181      REAL*8,save :: botCH4
    … …  
    10383      real*8,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
    10484      REAL*8,save,allocatable :: krpd(:,:,:,:),krate(:,:)
    105       integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
    106       integer,save :: prod(200,NC),perte(2,200,NC)
    107       character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
     85      integer,save :: reactif(5,nr_kim),nom_prod(nkim),nom_perte(nkim)
     86      integer,save :: prod(200,nkim),perte(2,200,nkim)
     87      character  fich*7,ficdec*15,curdec*15,name*10
    10888      real*8  ficalt,oldq,newq,zalt
    10989      logical okfic
    11090
     91!     Dummy parameters without microphysics
     92!      + Just here to keep the whole stuff running without modif C sources
     93!     ---------------------------------------------------------------------
     94
     95      INTEGER  :: utilaer(16)
     96      INTEGER  :: aerprod = 0
     97      INTEGER  :: htoh2   = 0
    11198
    11299!-----------------------------------------------------------------------
    … …  
    127114! ************************************
    128115
    129         allocate(kedd(NLEV))
    130 
    131         allocate(krpd(15,ND+1,NLRT,nbp_lat),krate(NLEV,NR),md(NLEV,NC))
    132 
    133 ! Verification du nombre de composes: coherence common_mod et nqmax-nmicro
    134 ! ----------------------------------
    135 
    136         if (ny.ne.NC) then
    137            print*,'PROBLEME de coherence dans le nombre de composes:',ny,NC
    138            STOP
    139         endif
     116        allocate(kedd(nlaykim_tot))
     117
     118        allocate(krpd(15,nd_kim+1,nlrt_kim,nbp_lat),krate(nlaykim_tot,nr_kim),md(nlaykim_tot,nkim))
     119
    140120
    141121! calcul de temp_c, densites et press_c en moyenne planetaire :
    … …  
    157137        rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
    158138
    159 ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
    160        do l=klev+2,NLEV
     139! au-dessus du GCM, dr regulier et rinter(nlaykim_tot)=1290+2575 km.
     140       do l=klev+2,nlaykim_tot
    161141         rinter(l) = rinter(klev+1) &
    162                + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
     142               + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
    163143         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
    164144       enddo
    165        rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
     145       rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-1))/2.
    166146! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire
    167147
    … …  
    169149! remplissage r1d,t1d,ct1d,p1d
    170150        open(11,file=TRIM(datadir)//'/tcp.ver',status='old')
    171         read(11,*) dummy
     151        read(11,*)
    172152        do i=1,131
    173153          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
    … …  
    176156        close(11)
    177157
    178 ! extension pour klev+1 a NLEV avec tcp.ver
     158! extension pour klev+1 a nlaykim_tot avec tcp.ver
    179159! la jonction klev/klev+1 est brutale... Tant pis ?
    180160        ialt=1
    181         do l=klev+1,NLEV
     161        do l=klev+1,nlaykim_tot
    182162           zalt=rmil(l)-rad/1000.
    183163           do i=ialt,130
    … …  
    200180        call comp(nomqy_c,nb,temp_c,mass,md)
    201181        print*,'           Mass'
    202         do ic=1,NC
     182        do ic=1,nkim
    203183          print*,nomqy_c(ic),mass(ic)
    204184!         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
    … …  
    214194                    nom_perte,nom_prod,perte,prod)
    215195!        print*,'nom_prod, nom_perte:'
    216 !        do ic=1,NC
     196!        do ic=1,nkim
    217197!          print*,nom_prod(ic),nom_perte(ic)
    218198!        enddo
    219199!        print*,'premieres prod, perte(1:reaction,2:compagnon):'
    220 !        do ic=1,NC
     200!        do ic=1,nkim
    221201!          print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
    222202!        enddo
    … …  
    224204!       l = klev-3
    225205!       print*,'krate a p=',press_c(l),' reactifs et produits:'
    226 !       do ic=1,ND+1
     206!       do ic=1,nd_kim+1
    227207!         print*,ic,krpd(7,ic,l,4)*nb(l),"  ",
    228208!    .     nomqy_c(reactif(1,ic)+1),
    … …  
    230210!    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
    231211!       enddo
    232 !       do ic=ND+2,NR
     212!       do ic=nd_kim+2,nr_kim
    233213!         print*,ic,krate(l,ic),"  ",
    234214!    .     nomqy_c(reactif(1,ic)+1),
    … …  
    245225      kedd(:) = 5e2 ! valeur mise par defaut 
    246226               ! UNITE ? doit etre ok pour gptitan
    247       do l=1,NLEV
     227      do l=1,nlaykim_tot
    248228       zalt=rmil(l)-rad/1000.  ! z en km
    249229       if     ((zalt.ge.250.).and.(zalt.le.400.)) then
    … …  
    280260!                * Permet de faire le calcul une seule fois par lat
    281261!
    282       DO j=1,nlon
     262      DO j=1,ngrid
    283263     
    284264      if (j.eq.1) then
    … …  
    305285       rinter(klev+1)=(rad+czlev(j,klev+1))/1000.
    306286
    307 ! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
    308        do l=klev+2,NLEV
     287! au-dessus du GCM, dr regulier et rinter(nlaykim_tot)=1290+2575 km.
     288       do l=klev+2,nlaykim_tot
    309289         rinter(l) = rinter(klev+1)  &
    310                + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
     290               + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
    311291         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
    312292       enddo
    313        rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
     293       rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-1))/2.
    314294
    315295!-----------------------------------------------------------------------
    … …  
    326306               nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l))
    327307       ENDDO
    328 ! extension pour klev+1 a NLEV avec tcp.ver
     308! extension pour klev+1 a nlaykim_tot avec tcp.ver
    329309       ialt=1
    330        do l=klev+1,NLEV
     310       do l=klev+1,nlaykim_tot
    331311           zalt=rmil(l)-rad/1000.
    332312           do i=ialt,130
    … …  
    371351      endif
    372352
    373       do l=1,NLEV
     353      do l=1,nlaykim_tot
    374354
    375355! INTERPOL EN ALT POUR k (krpd tous les deux km)
    … …  
    377357        factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
    378358
    379         do i=1,ND+1
     359        do i=1,nd_kim+1
    380360          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt) &
    381361                   + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
    … …  
    385365                   + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
    386366
    387                     ! ND+1 correspond a la dissociation de N2 par les GCR
     367                    ! nd_kim+1 correspond a la dissociation de N2 par les GCR
    388368          if ( factdec.lt.0. ) then
    389369        krate(l,i) = krpddecm1 * abs(factdec)  &
    … …  
    403383!   ------------
    404384
    405        do ic=1,NC
     385       do ic=1,nkim
    406386        do l=1,klev
    407387          cqy(l,ic)  = qy_c(j,l,ic)
    … …  
    410390       enddo
    411391
    412 ! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
     392! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a nlaykim_tot
    413393
    414394      WRITE(str2,'(i2.2)') klat(j)
    … …  
    423403! on lit la colonne t-1 au lieu de la colonne t
    424404! (cas d une bande de latitude partagee par 2 procs)
    425        do ic=1,NC
     405       do ic=1,nkim
    426406        read(11,'(A10)') name
    427407        if (name.ne.nomqy_c(ic)) then
    … …  
    431411        endif
    432412        if (ficdec.eq.curdec) then
    433           do l=klev+1,NLEV
     413          do l=klev+1,nlaykim_tot
    434414            read(11,*) ficalt,cqy(l,ic),newq
    435415          enddo
    436416        else
    437           do l=klev+1,NLEV
     417          do l=klev+1,nlaykim_tot
    438418            read(11,*) ficalt,oldq,cqy(l,ic)
    439419          enddo
    … …  
    442422       close(11)
    443423      else       ! le fichier n'est pas la
    444        do ic=1,NC
    445         do l=klev+1,NLEV
     424       do ic=1,nkim
     425        do l=klev+1,nlaykim_tot
    446426          cqy(l,ic)=cqy(klev,ic)
    447427        enddo
    … …  
    466446!   ---------------------
    467447     
    468        do ic=1,NC
     448       do ic=1,nkim
    469449        do l=1,klev
    470450          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
    … …  
    475455!-----------------------------------------------------------------------
    476456!
    477 !   sauvegarde compo sur NLEV
     457!   sauvegarde compo sur nlaykim_tot
    478458!   -------------------------
    479459
    … …  
    483463! premiere ligne=declin
    484464      write(11,'(E15.5)') declin_c
    485       do ic=1,NC
     465      do ic=1,nkim
    486466       write(11,'(A10)') nomqy_c(ic)
    487        do l=klev+1,NLEV
     467       do l=klev+1,nlaykim_tot
    488468        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-rad/1000.,cqy0(l,ic),cqy(l,ic)
    489469       enddo

  • trunk/LMDZ.TITAN/libf/phytitan/chem_settings.F90

    r1895 r1903  
    8383    ! ( H=1, H2=2 ..., C4N2=44) -> cf comchem_h
    8484   
    85     DO iq=2,44
     85    DO iq=2,nkim
    8686      CALL get_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:),found,indextime)
    8787      IF (.NOT.found) THEN

  • trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90

    r1899 r1903  
    2121IMPLICIT NONE 
    2222
     23   !! Hard-coded number of chemical species for Titan chemistry
     24   INTEGER, PARAMETER :: nkim = 44
     25
    2326   !! Hard-coded chemical species for Titan chemistry
    24    CHARACTER(len=10), DIMENSION(44), PARAMETER  :: cnames = &
     27   CHARACTER(len=10), DIMENSION(nkim), PARAMETER  :: cnames = &
    2528     (/"H         ", "H2        ", "CH        ", "CH2s      ", "CH2       ", "CH3       ", &
    2629       "CH4       ", "C2        ", "C2H       ", "C2H2      ", "C2H3      ", "C2H4      ", &
    … …  
    3134       "H2CN      ", "CHCN      ", "CH2CN     ", "CH3CN     ", "C3N       ", "HC3N      ", &
    3235       "NCCN      ", "C4N2      "/)
     36   !! Hard-coded chemical species for Titan chemistry + "HV" specie for the photochem module.
     37   CHARACTER(len=10), DIMENSION(nkim+1), PARAMETER  :: nomqy_c = &
     38     (/"H         ", "H2        ", "CH        ", "CH2s      ", "CH2       ", "CH3       ", &
     39       "CH4       ", "C2        ", "C2H       ", "C2H2      ", "C2H3      ", "C2H4      ", &
     40       "C2H5      ", "C2H6      ", "C3H3      ", "C3H5      ", "C3H6      ", "C3H7      ", &
     41       "C4H       ", "C4H3      ", "C4H4      ", "C4H2s     ", "CH2CCH2   ", "CH3CCH    ", &
     42       "C3H8      ", "C4H2      ", "C4H6      ", "C4H10     ", "AC6H6     ", "C3H2      ", &
     43       "C4H5      ", "AC6H5     ", "N2        ", "N4S       ", "CN        ", "HCN       ", &
     44       "H2CN      ", "CHCN      ", "CH2CN     ", "CH3CN     ", "C3N       ", "HC3N      ", &
     45       "NCCN      ", "C4N2      ", "HV        "/)
    3346   !! Hard-coded chemical species molar mass (g.mol-1), shares the same indexing than cnames.
    34    REAL, DIMENSION(44), PARAMETER               :: cmmol = (/ &
     47   REAL, DIMENSION(nkim), PARAMETER               :: cmmol = (/ &
    3548       1.01   , 2.0158, 13.02, 14.03, 14.03, 15.03, 16.04  , 24.02, 25.03, 26.04  , 27.05  , &
    3649       28.05  , 29.06 , 30.07, 39.06, 41.07, 42.08, 43.09  , 49.05, 51.07, 52.08  , 50.06  , &
    … …  
    3851       14.01  , 26.02 , 27.04, 28.05, 39.05, 40.04, 41.05  , 50.04, 51.05, 52.04  , 76.1   /)
    3952   
    40    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     53   ! !!!!!!!!!!!!!!!!!!!!
    4154   !  Upper chemistry
    42    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     55   ! !!!!!!!!!!!!!!!!!!!!
     56   
     57   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     58   ! These parameters as well as nkim above, MUST match titan.h in chimtitan
     59   INTEGER, PARAMETER :: nd_kim   = 54     ! Number of photodissociations
     60   INTEGER, PARAMETER :: nr_kim   = 377    ! Number of reactions in chemistry scheme
     61   INTEGER, PARAMETER :: nlrt_kim = 650    ! For the UV rad. transf., 650 levels of 2 km
     62   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4363   
    4464   INTEGER, SAVE :: nlaykim_up   ! Number of upper atm. layers for chemistry from GCM top to 4.5E-5 Pa (1300km)
    … …  
    7191    IF (.NOT.allocated(preskim)) ALLOCATE(preskim(nlaykim_up))
    7292    IF (.NOT.allocated(zlay_kim)) ALLOCATE(zlay_kim(ngrid,nlaykim_tot))
    73     IF (.NOT.allocated(ykim_up)) ALLOCATE(ykim_up(44,ngrid,nlaykim_up))
     93    IF (.NOT.allocated(ykim_up)) ALLOCATE(ykim_up(nkim,ngrid,nlaykim_up))
    7494 
    7595  END SUBROUTINE ini_comchem_h

  • trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90

    r1672 r1903  
    1       SUBROUTINE inicondens(ny,press,temp,nomy,yc)
     1      SUBROUTINE inicondens(press,temp,yc)
    22
    33!==================================================================================
    44!        Purpose
    55!        -------
    6 !        Initialisation des profils de saturation des traceurs
     6!        Initialisation des profils de saturation des traceurs chimiques
    77!
    88!        Authors
    … …  
    1717!   -------------
    1818
     19      use comchem_h, only: nkim, cnames
    1920      use dimphy
    2021      use mod_grid_phy_lmdz, only: nbp_lev
    … …  
    2324!    Arguments :
    2425!    -----------
    25       INTEGER, INTENT(IN)     :: ny
    26       CHARACTER*10,INTENT(IN) :: nomy(ny+1)
    27       REAL, INTENT(IN)        :: press(nbp_lev),temp(nbp_lev)  ! pressure in mbar !
    28       REAL, INTENT(OUT)       :: yc(nbp_lev,ny)
     26      REAL, INTENT(IN)        :: press(nbp_lev),temp(nbp_lev)  ! Pressure (mbar)
     27      REAL, INTENT(OUT)       :: yc(nbp_lev,nkim) ! Saturation profiles (mol/mol)
    2928     
    3029!    Local variables :
    … …  
    3332      REAL                    ::  sy,x
    3433       
    35       do ic=1,ny
    36        print*, 'traceur CH(', ic, ')=', nomy(ic),'------------'
     34      do ic=1,nkim
     35       print*, 'traceur CH(', ic, ')=', trim(cnames(ic)),'------------'
    3736           do l=1,nbp_lev
    3837
    … …  
    4039              yc(l,ic)=1.
    4140
    42               if(nomy(ic).eq."CH4") then
     41              if(trim(cnames(ic)).eq."CH4") then
    4342                 if (temp(l).lt.90.65) then
    4443                   yc(l,ic)= 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / TEMP(l) - &
    … …  
    5453              endif
    5554
    56               if(nomy(ic).eq."C2H2") then
     55              if(trim(cnames(ic)).eq."C2H2") then
    5756                 yc(l,ic)= 10.0**(6.09748e0-1644.1e0/TEMP(l)+7.42346e0     &
    5857                      * alog10(1.0e3/TEMP(l)) ) / PRESS(l)*1013.25e0/760.0
    5958              endif
    6059
    61               if(nomy(ic).eq."C2H4") then
     60              if(trim(cnames(ic)).eq."C2H4") then
    6261                 if (temp(l).lt.89.0) then
    6362                   yc(l,ic)= 10.0**(1.5477e0 + (1.0e0/TEMP(l) - 0.011e0)     &
    … …  
    7675              endif
    7776
    78               if(nomy(ic).eq."C2H6") then
     77              if(trim(cnames(ic)).eq."C2H6") then
    7978                 if (temp(l).lt.90.) then
    8079                   yc(l,ic)= 10.0**(10.01e0-1085.0e0/(TEMP(l)-0.561e0) )  &
    … …  
    8685              endif
    8786
    88               if((nomy(ic).eq."CH3CCH").or.(nomy(ic).eq."CH2CCH2")) then
     87              if((trim(cnames(ic)).eq."CH3CCH").or.(trim(cnames(ic)).eq."CH2CCH2")) then
    8988                 yc(l,ic)= 10.0**(2.8808e0 - 4.5e0*(249.9e0 - TEMP(l))   &
    9089                      /(1.15e0*TEMP(l) - 37.485e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
    9190              endif
    9291
    93               if(nomy(ic).eq."C3H6")  then
     92              if(trim(cnames(ic)).eq."C3H6")  then
    9493                 yc(l,ic)= 10.0**(7.4463e0 - 1028.5654e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    9594              endif
    9695
    97               if(nomy(ic).eq."C3H8")  then
     96              if(trim(cnames(ic)).eq."C3H8")  then
    9897                 yc(l,ic)= 10.0**(7.217e0 - 994.30251e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    9998              endif
    10099
    101               if((nomy(ic).eq."C4H2").or.(nomy(ic).eq."C4H2s")) then
     100              if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
    102101                 yc(l,ic)= 10.0**(96.26781e0 - 4651.872e0/TEMP(l) - 31.68595e0 &
    103102                      *alog10(TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
    104103              endif
    105104
    106               if(nomy(ic).eq."C4H4")  then
     105              if(trim(cnames(ic)).eq."C4H4")  then
    107106                 yc(l,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(TEMP(l)-43.15e0) ) / PRESS(l)
    108107              endif
    109108
    110               if(nomy(ic).eq."C4H6")  then
     109              if(trim(cnames(ic)).eq."C4H6")  then
    111110                 yc(l,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - TEMP(l)) &
    112111                      /(1.15e0*TEMP(l) - 39.345e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
    113112              endif
    114113
    115               if(nomy(ic).eq."C4H10")  then
     114              if(trim(cnames(ic)).eq."C4H10")  then
    116115                 yc(l,ic)= 10.0**(8.446e0 - 1461.2e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    117116              endif
    118117
    119               if(nomy(ic).eq."C6H2")  then
     118              if(trim(cnames(ic)).eq."C6H2")  then
    120119                 yc(l,ic)= 10.0**(4.666e0 - 4956e0/TEMP(l) + 25.845e0 * &
    121120                      alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
    122121              endif
    123122
    124               if(nomy(ic).eq."C8H2")  then
     123              if(trim(cnames(ic)).eq."C8H2")  then
    125124                 yc(l,ic)= 10.0**(3.95e0 - 6613e0/TEMP(l) + 35.055e0 * &
    126125                      alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
    127126              endif
    128127
    129               if(nomy(ic).eq."AC6H6")  then
     128              if(trim(cnames(ic)).eq."AC6H6")  then
    130129                 x = 1.0e0 - TEMP(l) / 562.2e0
    131130                 yc(l,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x &
    … …  
    133132              endif
    134133
    135               if(nomy(ic).eq."HCN")  then
     134              if(trim(cnames(ic)).eq."HCN")  then
    136135                 yc(l,ic)= 10.0**(8.6165e0 - 1516.5e0/(TEMP(l) - 26.2e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
    137136              endif
    138137
    139               if(nomy(ic).eq."CH3CN")  then
     138              if(trim(cnames(ic)).eq."CH3CN")  then
    140139                 yc(l,ic)= 10.0**(8.458e0 - 1911.7e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    141140              endif
    142141
    143               if(nomy(ic).eq."C2H3CN")  then
     142              if(trim(cnames(ic)).eq."C2H3CN")  then
    144143                 yc(l,ic)= 10.0**(9.3051e0 - 2782.21/(TEMP(l) - 51.15e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
    145144              endif
    146145
    147               if(nomy(ic).eq."NCCN")  then
     146              if(trim(cnames(ic)).eq."NCCN")  then
    148147                 yc(l,ic)=  10.0**(7.454e0 - 1832e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    149148              endif
    150149
    151               if(nomy(ic).eq."HC3N")  then
     150              if(trim(cnames(ic)).eq."HC3N")  then
    152151                 yc(l,ic)= 10.0**(7.7446e0 - 1453.5609e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    153152              endif
    154153
    155               if(nomy(ic).eq."C4N2")  then
     154              if(trim(cnames(ic)).eq."C4N2")  then
    156155                 yc(l,ic)= 10.0**(8.269e0 - 2155.0e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
    157156              endif

  • trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90

    r1894 r1903  
    143143  use iostart, only : open_restartphy, close_restartphy, &
    144144                      put_var, put_field
    145   use comchem_h, only: cnames, ykim_up
     145  use comchem_h, only: nkim, cnames, ykim_up
    146146  use tracer_h, only: noms
    147147  use callkeys_mod, only: callchim
    … …  
    196196  ! Upper chemistry
    197197  if (callchim) then
    198     do iq=1,44
     198    do iq=1,nkim
    199199      call put_field(trim(cnames(iq))//"_up",trim(cnames(iq))//" in upper atmosphere",ykim_up(iq,:,:))
    200200    enddo

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r1902 r1903  
    1919      use radcommon_h, only: sigma, glat, grav, BWNV
    2020      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
     21      use comchem_h, only: nkim
    2122      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
    2223      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
    … …  
    363364      character(len=10) :: tmp2
    364365
    365 ! Local variables for Titan chemistry and microphysics (JVO 2017)
    366 ! ----------------------------------------------------------------------------
    367      
    368       character*10,dimension(:),allocatable,save :: nomqy
    369 
    370 !$OMP THREADPRIVATE(nomqy)
    371 
     366! Local variables for Titan chemistry and microphysics
     367! ----------------------------------------------------
     368 
    372369      real ctimestep ! Chemistry timestep (s)
    373370 
    374       real temp_eq(nlayer), press_eq(nlayer) ! Moyennes planétaires
    375      
    376       real , allocatable, dimension(:,:,:),save :: ychim
    377 
    378       ! 2D vmr tendencies ( chemistry and condensation ) for chem. tracers
    379       real,dimension(:,:,:),allocatable,save :: dycchi ! Saved since chemistry is not called every step
    380       real dyccond(ngrid,nlayer,nq)       
    381          
    382       real,dimension(:,:),allocatable,save   :: qysat
    383      
    384 !$OMP THREADPRIVATE(ychim,dycchi,qysat)
    385 
    386       real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank
     371      ! Chemical tracers in molar fraction
     372      real, dimension(:,:,:), allocatable, save   :: ychim ! (mol/mol)
     373!$OMP THREADPRIVATE(ychim)
     374
     375      ! Molar fraction tendencies ( chemistry and condensation ) for tracers (mol/mol/s)
     376      real, dimension(ngrid,nlayer,nq)            :: dyccond ! All tracers, we want to use indx on it.
     377      real, dimension(:,:,:), allocatable, save   :: dycchi  ! Only for chem tracers. Saved since chemistry is not called every step.
     378!$OMP THREADPRIVATE(dycchi)
     379
     380      ! Saturation profiles
     381      real, dimension(:,:), allocatable, save     :: qysat ! (mol/mol)
     382!$OMP THREADPRIVATE(qysat)
     383      real temp_eq(nlayer), press_eq(nlayer) ! Planetary averages for the init. of saturation profiles.
     384
     385      ! Surface methane tank
     386      real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank (m)
    387387!$OMP THREADPRIVATE(tankCH4)
    388388
    … …  
    429429
    430430        ! Initialisation of nmicro as well as tracers names, indexes ...
    431         if (ngrid.ne.1) then
    432            call initracer2(nq,nametrac) ! Already done in rcm1d
     431        if (ngrid.ne.1) then ! Already done in rcm1d
     432           call initracer2(nq,nametrac)
    433433        endif
    434434
    … …  
    463463        ALLOCATE(zdtlw(ngrid,nlayer))
    464464        ALLOCATE(zdtsw(ngrid,nlayer))
    465         ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
    466         ALLOCATE(qysat(nlayer,nq-nmicro))
    467         ALLOCATE(nomqy(nq-nmicro+1)) ! +1 because of hv
    468        
     465 
    469466        ! This is defined in comsaison_h
    470467        ALLOCATE(mu0(ngrid))
    … …  
    506503!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    507504
    508          if ( callchim .and. (nq.gt.nmicro) ) then
    509 
    510             allocate(ychim(ngrid,nlayer,nq-nmicro))
    511 
     505         if ( callchim ) then
     506
     507            allocate(ychim(ngrid,nlayer,nkim))
     508            allocate(dycchi(ngrid,nlayer,nkim)) ! only for chemical tracers
     509            allocate(qysat(nlayer,nkim))
     510           
     511            ! Chemistry timestep
    512512            ctimestep = ptimestep*REAL(ichim)
    513513
    514             do iq=nmicro+1,nq
    515                nomqy(iq-nmicro) = nametrac(iq)
    516             enddo
    517 
    518             nomqy(nq-nmicro+1) = "HV"
    519            
    520             ! qysat is taken at the equator ( small variations of t,p)
     514            ! qysat is taken at the equator ( small variations of t,p )
    521515            temp_eq(:)  = tmoy(:)
    522             press_eq(:) = playmoy(:)/100. ! en mbar
    523            
    524             call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
     516            press_eq(:) = playmoy(:)/100. ! in mbar
     517           
     518            call inicondens(press_eq,temp_eq,qysat)
    525519         
    526520            zdqchi(:,:,:)  = 0.0
    … …  
    10841078         if (callchim) then
    10851079
    1086             if (nq.gt.nmicro) then
    1087                do iq = nmicro+1,nq
    1088                   ychim(:,:,iq-nmicro) = pq(:,:,iq) * rat_mmol(iq) ! convert to molar fraction
    1089                enddo
    1090             endif
     1080            do iq = 1,nkim
     1081               ychim(:,:,iq) = pq(:,:,iq+nmicro) * rat_mmol(iq+nmicro) ! convert to molar fraction
     1082            enddo
    10911083
    10921084            ! Condensation tendency after the transport
    1093             do iq=nmicro+1,nq
     1085            do iq=1,nkim
    10941086               do l=1,nlayer
    10951087                  do ig=1,ngrid
    1096                      if ( ychim(ig,l,iq-nmicro).gt.qysat(l,iq-nmicro) ) then
    1097                         dyccond(ig,l,iq)= ( -ychim(ig,l,iq-nmicro)+qysat(l,iq-nmicro) ) / ptimestep
     1088                     if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
     1089                        dyccond(ig,l,iq+nmicro)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
    10981090                     endif
    10991091                  enddo
    … …  
    11081100
    11091101               if (ngrid.eq.1) then ! We obviously don't have access to (and don't need) zonal means in 1D
    1110                  call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
    1111                            pt,pplay,pplev,zzlay,zzlev,dycchi,nlayer+70)
     1102                 call calchim(ngrid,ychim,declin,zls,ctimestep, &
     1103                           pt,pplay,pplev,zzlay,zzlev,dycchi)
    11121104               else
    1113                  call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
    1114                            ztfibar,zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi,nlayer+70)
    1115                  ! JVO 2017 : NLEV = nlayer+70, en accord avec le C. Quid si nlay=/ 55 ?
     1105                 call calchim(ngrid,ychim,declin,zls,ctimestep, &
     1106                           ztfibar,zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
    11161107               endif
    11171108
    … …  
    11221113            ! ( GCM-friendly tracers and tendencies -> format for photochem and microphys )
    11231114
    1124             if (nq.gt.nmicro) then
    1125                ! We convert tendencies to mass mixing ratio
    1126                do iq=nmicro+1,nq
    1127                   zdqchi(:,:,iq)  = dycchi(:,:,iq-nmicro)  / rat_mmol(iq)
    1128                   zdqcond(:,:,iq) = dyccond(:,:,iq) / rat_mmol(iq)
    1129                enddo
    1130 
    1131                pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
    1132                     zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
     1115            ! We convert tendencies to mass mixing ratio
     1116            do iq=1,nkim
     1117               zdqchi(:,:,iq+nmicro)  = dycchi(:,:,iq)  / rat_mmol(iq+nmicro)
     1118               zdqcond(:,:,iq+nmicro) = dyccond(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
     1119            enddo
     1120
     1121            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
     1122                 zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
    11331123               
    1134             endif
    11351124           
    11361125         endif ! end of 'callchim'

  • trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90

    r1894 r1903  
    6868    USE callkeys_mod
    6969    USE comcstfi_mod, only: mugaz
    70     USE comchem_h, only: cnames, cmmol
     70    USE comchem_h, only: nkim, cnames, cmmol
    7171    IMPLICIT NONE
    7272
    … …  
    146146
    147147    IF (callchim) THEN
    148       IF (nchimi < SIZE(cnames)) THEN
    149         WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
     148      IF (nchimi .NE. nkim) THEN
     149        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",nkim," expected)"
    150150        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
    151151      ENDIF
    152152      IF (.NOT.ALLOCATED(chimi_indx)) ALLOCATE(chimi_indx(nchimi))
    153153      n = 0 ! counter on chimi_indx
    154       DO j=1,SIZE(cnames)
     154      DO j=1,nkim
    155155        found = .false.
    156156        DO i=1,nq
    … …  
    168168        ENDIF
    169169      ENDDO
    170       IF (n < SIZE(cnames)) THEN
    171         WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
     170      IF (n .NE. nkim) THEN
     171        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",nkim," expected)"
    172172        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
    173173      ENDIF
Note: See TracChangeset for help on using the changeset viewer.