Ignore:
Timestamp:
Apr 13, 2012, 8:40:01 AM (13 years ago)
Author:
emillour
Message:

Mars GCM:

Update of the chemistry package:

  • calchim.F90 :
    • upgraded from calchim.F
    • can run with or without CH4
    • fixed the mass conservation scheme
  • photochemistry.F : -chemistry timestep is now independant from physics timestep -can run with or without a CH4 tracer
    • removed initial tests on species, since these are already done in calchim
  • removed inichim_readcallphys.F (not used any more)
  • concentrations.F: adaptated to better handle indexes and molar masses of

tracers. "ncomp" (in chimiedata.h) no longer needed.

  • inichim_newstart.F90 : rewritten and cleaned (won't be compatible with

old start files where tracer names are q01,...).
Now handles hybrid levels and automaticaly adapts
depending on which tracers are available.

  • newstart.F : adapted to followup changes in inchim_newstart.F90 and some

cleanup around the initialization of tracer names and surface
values.

FL+EM

Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
1 deleted
2 edited
2 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim.F90

    r616 r618  
    1       subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
    2      $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
    3      $                   dqscloud,tauref,co2ice,
    4      $                   pu,pdu,pv,pdv,surfdust,surfice)
    5 c
     1      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,         &
     2                         zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,    &
     3                         dqscloud,tauref,co2ice,                            &
     4                         pu,pdu,pv,pdv,surfdust,surfice)
     5
    66      implicit none
    7 c
    8 c=======================================================================
    9 c
    10 c   subject:
    11 c   --------
    12 c
    13 c  Prepare the call for the photochemical module, and send back the
    14 c  tendencies from photochemistry in the chemical species mass mixing ratios
    15 c
    16 c   Author:   Sebastien Lebonnois (08/11/2002)
    17 c   -------
    18 c    update 12/06/2003 for water ice clouds and compatibility with dust
    19 c    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
    20 c    update 03/05/2005 cosmetic changes (Franck Lefevre)
    21 c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
    22 c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
    23 c    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
    24 c
    25 c   Arguments:
    26 c   ----------
    27 c
    28 c  Input:
    29 c
    30 c    ptimestep                  timestep (s)
    31 c    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
    32 c    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
    33 c    pt(ngridmx,nlayermx)       Temperature (K)
    34 c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
    35 c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
    36 c    pdu(ngridmx,nlayermx)      u component tendency (K)
    37 c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
    38 c    pdv(ngridmx,nlayermx)      v component tendency (K)
    39 c    dist_sol                   distance of the sun (AU)
    40 c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
    41 c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
    42 c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
    43 c    tauref(ngridmx)            Optical depth at 7 hPa
    44 c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
    45 c    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
    46 c    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
    47 c
    48 c  Output:
    49 c
    50 c    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
    51 c    dqschim(ngridmx,nqmx)         ! tendencies on qsurf
    52 c
    53 c=======================================================================
    54 
    55 c    Declarations :
    56 c    --------------
     7
     8!=======================================================================
     9!
     10!   subject:
     11!   --------
     12!
     13!  Prepare the call for the photochemical module, and send back the
     14!  tendencies from photochemistry in the chemical species mass mixing ratios
     15!
     16!   Author:   Sebastien Lebonnois (08/11/2002)
     17!   -------
     18!    update 12/06/2003 for water ice clouds and compatibility with dust
     19!    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
     20!    update 03/05/2005 cosmetic changes (Franck Lefevre)
     21!    update sept. 2008 identify tracers by their names (Ehouarn Millour)
     22!    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
     23!    update 16/03/2012 optimization (Franck Lefevre)
     24!
     25!   Arguments:
     26!   ----------
     27!
     28!  Input:
     29!
     30!    ptimestep                  timestep (s)
     31!    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
     32!    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
     33!    pt(ngridmx,nlayermx)       Temperature (K)
     34!    pdt(ngridmx,nlayermx)      Temperature tendency (K)
     35!    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
     36!    pdu(ngridmx,nlayermx)      u component tendency (K)
     37!    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
     38!    pdv(ngridmx,nlayermx)      v component tendency (K)
     39!    dist_sol                   distance of the sun (AU)
     40!    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
     41!    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
     42!    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
     43!    tauref(ngridmx)            Optical depth at 7 hPa
     44!    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
     45!    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
     46!    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
     47!
     48!  Output:
     49!
     50!    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
     51!    dqschim(ngridmx,nqmx)         ! tendencies on qsurf
     52!
     53!=======================================================================
    5754
    5855#include "dimensions.h"
    … …  
    6461#include "conc.h"
    6562
    66 c Arguments :
    67 c -----------
    68 
    69 c   inputs:
    70 c   -------
    71 
    72       real    ptimestep
    73       real    pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
    74       real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
    75       real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
    76       real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
    77       real    pt(ngridmx,nlayermx)       ! temperature
    78       real    pdt(ngridmx,nlayermx)      ! temperature tendency
    79       real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
    80       real    pdu(ngridmx,nlayermx)      ! u component tendency
    81       real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
    82       real    pdv(ngridmx,nlayermx)      ! v component tendency
    83       real    dist_sol                   ! distance of the sun (AU)
    84       real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
    85       real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
    86       real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
    87       real    zday                       ! date (time since Ls=0, in martian days)
    88       real    tauref(ngridmx)            ! optical depth at 7 hPa
    89       real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
    90       real    surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
    91       real    surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
    92 
    93 c   outputs:
    94 c   --------
    95 
    96       real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
    97       real dqschim(ngridmx,nqmx)         ! tendencies on qsurf
    98       real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
    99       real dqscloud(ngridmx,nqmx)        ! tendencies on qsurf
    100 
    101 c Local variables :
    102 c -----------------
    103 
    104       character*20 str20
    105       integer  ig,l,i,iq
    106       integer  foundswitch, lswitch
    107       integer  ig_vl1
    108       real latvl1, lonvl1
    109       real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
    110                                      ! new mole fraction after
    111       real zt(ngridmx,nlayermx)      ! temperature
    112       real zu(ngridmx,nlayermx)      ! u component of the wind
    113       real zv(ngridmx,nlayermx)      ! v component of the wind
    114       real taucol                    ! optical depth at 7 hPa
    115 
    116       logical depos                  ! switch for dry deposition
    117       parameter (depos=.false.)
    118 c
    119 c  for each column of atmosphere:
    120 c
    121       real zpress(nlayermx)       !  Pressure (mbar)
    122       real zdens(nlayermx)        !  Density  (cm-3)
    123       real ztemp(nlayermx)        !  Temperature (K)
    124       real zlocal(nlayermx)       !  Altitude (km)
    125       real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
    126       real szacol                 !  Solar zenith angle
    127       real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
    128       real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
    129       real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
    130 c
    131 c for output:
    132 c
    133       real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)
    134 
    135       logical output              ! to issue calls to writediagfi and stats
    136       parameter (output=.true.)   ! see at end of routine
    137 
    138       logical,save :: firstcall=.true.
    139       integer,save :: nbq       ! number of tracers used in the chemistry
    140       integer,save :: niq(nqmx) ! array storing the indexes of the tracers
    141 
    142 ! index of tracers:
    143 
    144       integer,save :: i_co2=0
    145       integer,save :: i_co=0
    146       integer,save :: i_o=0
    147       integer,save :: i_o1d=0
    148       integer,save :: i_o2=0
    149       integer,save :: i_o3=0
    150       integer,save :: i_h=0
    151       integer,save :: i_h2=0
    152       integer,save :: i_oh=0
    153       integer,save :: i_ho2=0
    154       integer,save :: i_h2o2=0
    155       integer,save :: i_ch4=0
    156       integer,save :: i_n2=0
    157       integer,save :: i_ar=0
    158       integer,save :: i_ice=0 ! water ice
    159       integer,save :: i_h2o=0 ! water vapour
    160 c
    161 c=======================================================================
    162 c     initialization of the chemistry (first call only)
    163 c=======================================================================
    164 c
     63!     input:
     64
     65      real :: ptimestep
     66      real :: pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
     67      real :: zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
     68      real :: pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
     69      real :: zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
     70      real :: pt(ngridmx,nlayermx)       ! temperature
     71      real :: pdt(ngridmx,nlayermx)      ! temperature tendency
     72      real :: pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
     73      real :: pdu(ngridmx,nlayermx)      ! u component tendency
     74      real :: pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
     75      real :: pdv(ngridmx,nlayermx)      ! v component tendency
     76      real :: dist_sol                   ! distance of the sun (AU)
     77      real :: mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
     78      real :: pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
     79      real :: pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
     80      real :: zday                       ! date (time since Ls=0, in martian days)
     81      real :: tauref(ngridmx)            ! optical depth at 7 hPa
     82      real :: co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
     83      real :: surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
     84      real :: surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)
     85
     86!     output:
     87
     88      real :: dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
     89      real :: dqschim(ngridmx,nqmx)         ! tendencies on qsurf
     90      real :: dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
     91      real :: dqscloud(ngridmx,nqmx)        ! tendencies on qsurf
     92
     93!     local variables:
     94
     95      integer,save :: nbq                   ! number of tracers used in the chemistry
     96      integer,save :: niq(nqmx)             ! array storing the indexes of the tracers
     97      integer :: major(nlayermx)            ! index of major species
     98      integer :: ig,l,i,iq,iqmax
     99      integer :: foundswitch, lswitch
     100
     101      integer,save :: i_co2  = 0
     102      integer,save :: i_co   = 0
     103      integer,save :: i_o    = 0
     104      integer,save :: i_o1d  = 0
     105      integer,save :: i_o2   = 0
     106      integer,save :: i_o3   = 0
     107      integer,save :: i_h    = 0
     108      integer,save :: i_h2   = 0
     109      integer,save :: i_oh   = 0
     110      integer,save :: i_ho2  = 0
     111      integer,save :: i_h2o2 = 0
     112      integer,save :: i_ch4  = 0
     113      integer,save :: i_n2   = 0
     114      integer,save :: i_h2o  = 0
     115
     116      integer :: ig_vl1
     117
     118      real    :: latvl1, lonvl1
     119      real    :: zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
     120                                           ! new mole fraction after
     121      real    :: zt(ngridmx,nlayermx)      ! temperature
     122      real    :: zu(ngridmx,nlayermx)      ! u component of the wind
     123      real    :: zv(ngridmx,nlayermx)      ! v component of the wind
     124      real    :: taucol                    ! optical depth at 7 hPa
     125
     126      logical,save :: firstcall = .true.
     127      logical,save :: depos = .false.      ! switch for dry deposition
     128
     129!     for each column of atmosphere:
     130
     131      real :: zpress(nlayermx)       !  Pressure (mbar)
     132      real :: zdens(nlayermx)        !  Density  (cm-3)
     133      real :: ztemp(nlayermx)        !  Temperature (K)
     134      real :: zlocal(nlayermx)       !  Altitude (km)
     135      real :: zycol(nlayermx,nqmx)   !  Composition (mole fractions)
     136      real :: szacol                 !  Solar zenith angle
     137      real :: surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
     138      real :: surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
     139      real :: jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
     140
     141!     for output:
     142
     143      logical :: output                 ! to issue calls to writediagfi and stats
     144      parameter (output = .true.)
     145      real :: jo3_3d(ngridmx,nlayermx)  ! Photodissociation rate O3->O1D (s-1)
     146
     147!=======================================================================
     148!     initialization of the chemistry (first call only)
     149!=======================================================================
     150
    165151      if (firstcall) then
    166 c
     152
    167153         if (photochem) then
    168154            print*,'calchim: Read photolysis lookup table'
    … …  
    171157
    172158         ! find index of chemical tracers to use
    173          nbq=0 ! to count number of tracers
    174          i_co2=igcm_co2
    175          if (i_co2.eq.0) then
    176            write(*,*) "calchim: Error; no CO2 tracer !!!"
    177            stop
    178          else
    179            nbq=nbq+1
    180            niq(nbq)=i_co2
    181          endif
    182          i_co=igcm_co
    183          if (i_co.eq.0) then
    184            write(*,*) "calchim: Error; no CO tracer !!!"
    185            stop
    186          else
    187            nbq=nbq+1
    188            niq(nbq)=i_co
    189          endif
    190          i_o=igcm_o
    191          if (i_o.eq.0) then
    192            write(*,*) "calchim: Error; no O tracer !!!"
    193            stop
    194          else
    195            nbq=nbq+1
    196            niq(nbq)=i_o
    197          endif
    198          i_o1d=igcm_o1d
    199          if (i_o1d.eq.0) then
    200            write(*,*) "calchim: Error; no O1D tracer !!!"
    201            stop
    202          else
    203            nbq=nbq+1
    204            niq(nbq)=i_o1d
    205          endif
    206          i_o2=igcm_o2
    207          if (i_o2.eq.0) then
    208            write(*,*) "calchim: Error; no O2 tracer !!!"
    209            stop
    210          else
    211            nbq=nbq+1
    212            niq(nbq)=i_o2
    213          endif
    214          i_o3=igcm_o3
    215          if (i_o3.eq.0) then
    216            write(*,*) "calchim: Error; no O3 tracer !!!"
    217            stop
    218          else
    219            nbq=nbq+1
    220            niq(nbq)=i_o3
    221          endif
    222          i_h=igcm_h
    223          if (i_h.eq.0) then
    224            write(*,*) "calchim: Error; no H tracer !!!"
    225            stop
    226          else
    227            nbq=nbq+1
    228            niq(nbq)=i_h
    229          endif
    230          i_h2=igcm_h2
    231          if (i_h2.eq.0) then
    232            write(*,*) "calchim: Error; no H2 tracer !!!"
    233            stop
    234          else
    235            nbq=nbq+1
    236            niq(nbq)=i_h2
    237          endif
    238          i_oh=igcm_oh
    239          if (i_oh.eq.0) then
    240            write(*,*) "calchim: Error; no OH tracer !!!"
    241            stop
    242          else
    243            nbq=nbq+1
    244            niq(nbq)=i_oh
    245          endif
    246          i_ho2=igcm_ho2
    247          if (i_ho2.eq.0) then
    248            write(*,*) "calchim: Error; no HO2 tracer !!!"
    249            stop
    250          else
    251            nbq=nbq+1
    252            niq(nbq)=i_ho2
    253          endif
    254          i_h2o2=igcm_h2o2
    255          if (i_h2o2.eq.0) then
    256            write(*,*) "calchim: Error; no H2O2 tracer !!!"
    257            stop
    258          else
    259            nbq=nbq+1
    260            niq(nbq)=i_h2o2
    261          endif
    262          i_ch4=igcm_ch4
    263          if (i_ch4.eq.0) then
    264            write(*,*) "calchim: Error; no CH4 tracer !!!"
    265            stop
    266          else
    267            nbq=nbq+1
    268            niq(nbq)=i_ch4
    269          endif
    270          i_n2=igcm_n2
    271          if (i_n2.eq.0) then
    272            write(*,*) "calchim: Error; no N2 tracer !!!"
    273            stop
    274          else
    275            nbq=nbq+1
    276            niq(nbq)=i_n2
    277          endif
    278          i_ar=igcm_ar
    279          if (i_ar.eq.0) then
    280            write(*,*) "calchim: Error; no AR tracer !!!"
    281            stop
    282          else
    283            nbq=nbq+1
    284            niq(nbq)=i_ar
    285          endif
    286          i_ice=igcm_h2o_ice
    287          if (i_ice.eq.0) then
    288            write(*,*) "calchim: Error; no water ice tracer !!!"
    289            stop
    290          else
    291            nbq=nbq+1
    292            niq(nbq)=i_ice
    293          endif
    294          i_h2o=igcm_h2o_vap
    295          if (i_h2o.eq.0) then
    296            write(*,*) "calchim: Error; no water vapor tracer !!!"
    297            stop
    298          else
    299            nbq=nbq+1
    300            niq(nbq)=i_h2o
    301          endif
     159         nbq = 0 ! to count number of tracers
     160         i_co2 = igcm_co2
     161         if (i_co2 == 0) then
     162            write(*,*) "calchim: Error; no CO2 tracer !!!"
     163            stop
     164         else
     165            nbq = nbq + 1
     166            niq(nbq) = i_co2
     167         end if
     168         i_co = igcm_co
     169         if (i_co == 0) then
     170            write(*,*) "calchim: Error; no CO tracer !!!"
     171            stop
     172         else
     173            nbq = nbq + 1
     174            niq(nbq) = i_co
     175         end if
     176         i_o = igcm_o
     177         if (i_o == 0) then
     178            write(*,*) "calchim: Error; no O tracer !!!"
     179            stop
     180         else
     181            nbq = nbq + 1
     182            niq(nbq) = i_o
     183         end if
     184         i_o1d = igcm_o1d
     185         if (i_o1d == 0) then
     186            write(*,*) "calchim: Error; no O1D tracer !!!"
     187            stop
     188         else
     189            nbq = nbq + 1
     190            niq(nbq) = i_o1d
     191         end if
     192         i_o2 = igcm_o2
     193         if (i_o2 == 0) then
     194            write(*,*) "calchim: Error; no O2 tracer !!!"
     195            stop
     196         else
     197            nbq = nbq + 1
     198            niq(nbq) = i_o2
     199         end if
     200         i_o3 = igcm_o3
     201         if (i_o3 == 0) then
     202            write(*,*) "calchim: Error; no O3 tracer !!!"
     203            stop
     204         else
     205            nbq = nbq + 1
     206            niq(nbq) = i_o3
     207         end if
     208         i_h = igcm_h
     209         if (i_h == 0) then
     210            write(*,*) "calchim: Error; no H tracer !!!"
     211            stop
     212         else
     213            nbq = nbq + 1
     214            niq(nbq) = i_h
     215         end if
     216         i_h2 = igcm_h2
     217         if (i_h2 == 0) then
     218            write(*,*) "calchim: Error; no H2 tracer !!!"
     219            stop
     220         else
     221            nbq = nbq + 1
     222            niq(nbq) = i_h2
     223         end if
     224         i_oh = igcm_oh
     225         if (i_oh == 0) then
     226            write(*,*) "calchim: Error; no OH tracer !!!"
     227            stop
     228         else
     229            nbq = nbq + 1
     230            niq(nbq) = i_oh
     231         end if
     232         i_ho2 = igcm_ho2
     233         if (i_ho2 == 0) then
     234            write(*,*) "calchim: Error; no HO2 tracer !!!"
     235            stop
     236         else
     237            nbq = nbq + 1
     238            niq(nbq) = i_ho2
     239         end if
     240         i_h2o2 = igcm_h2o2
     241         if (i_h2o2 == 0) then
     242            write(*,*) "calchim: Error; no H2O2 tracer !!!"
     243            stop
     244         else
     245            nbq = nbq + 1
     246            niq(nbq) = i_h2o2
     247         end if
     248         i_ch4 = igcm_ch4
     249         if (i_ch4 == 0) then
     250            write(*,*) "calchim: no CH4 tracer !!!"
     251            write(*,*) "CH4 will be ignored in the chemistry"
     252         else
     253            nbq = nbq + 1
     254            niq(nbq) = i_ch4
     255         end if
     256         i_n2 = igcm_n2
     257         if (i_n2 == 0) then
     258            write(*,*) "calchim: Error; no N2 tracer !!!"
     259            stop
     260         else
     261            nbq = nbq + 1
     262            niq(nbq) = i_n2
     263         end if
     264         i_h2o = igcm_h2o_vap
     265         if (i_h2o == 0) then
     266            write(*,*) "calchim: Error; no water vapor tracer !!!"
     267            stop
     268         else
     269            nbq = nbq + 1
     270            niq(nbq) = i_h2o
     271         end if
    302272
    303273         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
    304          write(*,*) '               i_co2  = ',i_co2
    305          write(*,*) '               i_co   = ',i_co
    306          write(*,*) '               i_o    = ',i_o
    307          write(*,*) '               i_o1d  = ',i_o1d
    308          write(*,*) '               i_o2   = ',i_o2
    309          write(*,*) '               i_o3   = ',i_o3
    310          write(*,*) '               i_h    = ',i_h
    311          write(*,*) '               i_h2   = ',i_h2
    312          write(*,*) '               i_oh   = ',i_oh
    313          write(*,*) '               i_ho2  = ',i_ho2
    314          write(*,*) '               i_h2o2 = ',i_h2o2
    315          write(*,*) '               i_ch4  = ',i_ch4
    316          write(*,*) '               i_n2   = ',i_n2
    317          write(*,*) '               i_ar   = ',i_ar
    318          write(*,*) '               i_ice  = ',i_ice
    319          write(*,*) '               i_h2o  = ',i_h2o
    320274         
    321275         firstcall = .false.
    322276      end if ! if (firstcall)
    323277
    324 ! Initialize output tendencies to zero (to handle case of tracers which
    325 ! are not used in the chemistry (e.g. dust))
    326 
     278! Initializations
     279
     280      zycol(:,:)    = 0.
    327281      dqchim(:,:,:) = 0
    328282      dqschim(:,:)  = 0
    329 c
     283
    330284!     latvl1= 22.27
    331285!     lonvl1= -47.94
    332 !     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
    333 !    $        int(1.5+(lonvl1+180)*iim/360.)
    334 c
    335 c=======================================================================
    336 c     loop over grid
    337 c=======================================================================
    338 c
     286!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +    &
     287!             int(1.5+(lonvl1+180)*iim/360.)
     288
     289!=======================================================================
     290!     loop over grid
     291!=======================================================================
     292
    339293      do ig = 1,ngridmx
    340 c
    341 c     local updates
    342 c
     294
    343295         foundswitch = 0
    344296         do l = 1,nlayermx
    345297            do i = 1,nbq
    346               iq = niq(i) ! get tracer index
    347               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
    348               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
     298               iq = niq(i) ! get tracer index
     299               zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
     300               zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
    349301            end do
    350             zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
    351             zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
    352             zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
    353 c
     302            zt(ig,l)  = pt(ig,l) + pdt(ig,l)*ptimestep
     303            zu(ig,l)  = pu(ig,l) + pdu(ig,l)*ptimestep
     304            zv(ig,l)  = pv(ig,l) + pdv(ig,l)*ptimestep
    354305            zpress(l) = pplay(ig,l)/100.
    355306            ztemp(l)  = zt(ig,l)
    356307            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
    357308            zlocal(l) = zzlay(ig,l)/1000.
    358 c
    359 c           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
    360 c
     309
     310!           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
     311
    361312            surfdust1d(l) = surfdust(ig,l)*1.e-2
    362313            surfice1d(l)  = surfice(ig,l)*1.e-2
    363 c
    364 c           search for switch index between regions
    365 c
     314
     315!           search for switch index between regions
     316
    366317            if (photochem .and. thermochem) then
    367                if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
     318               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
    368319                  lswitch = l
    369                   foundswitch=1
     320                  foundswitch = 1
    370321               end if
    371322            end if
    372             if ( .not. photochem) then
     323            if (.not. photochem) then
    373324               lswitch = 22
    374325            end if
    375             if (.not.  thermochem) then
    376                lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
     326            if (.not. thermochem) then
     327               lswitch = min(50,nlayermx+1)
    377328            end if
    378 c
     329
    379330         end do ! of do l=1,nlayermx
    380 c
     331
    381332         szacol = acos(mu0(ig))*180./pi
    382          taucol = tauref(ig)
    383 c
    384 c=======================================================================
    385 c     call chemical subroutine
    386 c=======================================================================
    387 c
    388 c        chemistry in lower atmosphere
    389 c
     333         taucol = tauref(ig)*(700./610.)  ! provisoire en attente de nouveau jmars
     334
     335!=======================================================================
     336!     call chemical subroutines
     337!=======================================================================
     338
     339!        chemistry in lower atmosphere
     340
    390341         if (photochem) then
    391             call photochemistry(lswitch,zycol,szacol,ptimestep,
    392      $                          zpress,ztemp,zdens,dist_sol,
    393      $                          surfdust1d,surfice1d,jo3,taucol)
    394          end if
    395 c
    396 c        chemistry in upper atmosphere
    397 c
     342            call photochemistry(lswitch,zycol,szacol,ptimestep,    &
     343                                zpress,ztemp,zdens,dist_sol,       &
     344                                surfdust1d,surfice1d,jo3,taucol)
     345
     346!        ozone photolysis, for output
     347
     348            do l = 1,nlayermx
     349               jo3_3d(ig,l) = jo3(l)
     350            end do
     351
     352!        condensation of h2o2
     353
     354            call perosat(ig,ptimestep,pplev,pplay,                 &
     355                         ztemp,zycol,dqcloud,dqscloud)
     356         end if
     357
     358!        chemistry in upper atmosphere
     359
    398360         if (thermochem) then
    399             call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
    400      $                       zlocal,szacol,ptimestep,zday)
    401          end if
    402 c
    403 c        dry deposition
    404 c
     361            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,  &
     362                             zlocal,szacol,ptimestep,zday)
     363         end if
     364
     365!        dry deposition
     366
    405367         if (depos) then
    406             call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,
    407      $                      zu, zv, zt, zycol, ptimestep, co2ice)
    408          end if
    409 c
    410 c=======================================================================
    411 c     tendencies
    412 c=======================================================================
    413 c
    414 c     must be 0. for water ice:
    415 c
    416          if (water) then
    417             do l = 1,nlayermx
    418                dqchim(ig,l,i_ice) = 0.
     368            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev,&
     369                            zu, zv, zt, zycol, ptimestep, co2ice)
     370         end if
     371
     372!=======================================================================
     373!     tendencies
     374!=======================================================================
     375
     376!     index of the most abundant species at each level
     377
     378         major(:) = maxloc(zycol, dim = 2)
     379
     380!     tendency for the most abundant species = - sum of others
     381
     382         do l = 1,nlayermx
     383            iqmax = major(l)
     384            do i = 1,nbq
     385               iq = niq(i) ! get tracer index
     386               if (iq /= iqmax) then
     387                  dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)  &
     388                                   - zq(ig,l,iq))/ptimestep
     389                  dqchim(ig,l,iqmax) = dqchim(ig,l,iqmax)              &
     390                                     - dqchim(ig,l,iq)
     391               end if
    419392            end do
    420          end if
    421 c
    422 c     tendency for CO2 = - sum of others for lower atmosphere
    423 c     tendency for O   = - sum of others for upper atmosphere
    424 c
    425          do l = 1,nlayermx
    426             if (l .lt. lswitch) then
    427                do i=1,nbq
    428                  iq=niq(i) ! get tracer index
    429                   if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then
    430                      dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)
    431      $                                - zq(ig,l,iq))/ptimestep
    432                   else if (iq.eq.i_co2) then
    433                      dqchim(ig,l,iq) = 0.
    434                   end if
    435                   dqschim(ig,iq) = 0.
    436                end do ! of do i=1,nbq
    437                
    438                do i=1,nbq
    439                  iq=niq(i) ! get tracer index
    440                   if (iq.ne.i_co2) then
    441                      dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2)
    442      $                                  - dqchim(ig,l,iq)
    443                   end if
    444                end do
    445              else if (l .ge. lswitch) then
    446                do i=1,nbq
    447                  iq=niq(i) ! get tracer index
    448                    if ((iq.ne.i_o).and.(iq.ne.i_ice)) then
    449                       dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)
    450      $                                  /mmean(ig,l)
    451      $                                 - zq(ig,l,iq))/ptimestep
    452                    else if (iq.eq.i_o) then
    453                       dqchim(ig,l,iq) = 0.
    454                    end if
    455                enddo
    456                
    457                do i=1,nbq
    458                  iq=niq(i) ! get tracer index
    459                    if (iq.ne.i_o) then
    460                       dqchim(ig,l,i_o) = dqchim(ig,l,i_o)
    461      $                                 - dqchim(ig,l,iq)
    462                    end if
    463                end do
    464              end if ! of if (l.lt.lswitch) else if (l.ge.lswitch)
    465           end do ! of do l = 1,nlayermx
    466 c
    467 c     condensation of h2o2
    468 c
    469           call perosat(ig,ptimestep,pplev,pplay,
    470      $                 ztemp,zycol,dqcloud,dqscloud)
    471 c
    472 c     density and j(o3->o1d), for outputs
    473 c
    474           do l = 1,nlayermx
    475              jo3_3d(ig,l) = jo3(l)
    476           end do
    477 c
    478 c=======================================================================
    479 c     end of loop over grid
    480 c=======================================================================
    481 c
     393         end do ! of do l = 1,nlayermx
     394
     395!=======================================================================
     396!     end of loop over grid
     397!=======================================================================
     398
    482399      end do ! of do ig=1,ngridmx
    483 c
    484 c=======================================================================
    485 c     write outputs
    486 c=======================================================================
    487 c
     400
     401!=======================================================================
     402!     write outputs
     403!=======================================================================
     404
    488405! value of parameter 'output' to trigger writting of outputs
    489406! is set above at the declaration of the variable.
    490407
    491       if (output) then
    492          if (ngridmx .gt. 1) then
    493             call writediagfi(ngridmx,'jo3','j o3->o1d',
    494      $                       's-1',3,jo3_3d(1,1))
     408      if (photochem .and. output) then
     409         if (ngridmx > 1) then
     410            call writediagfi(ngridmx,'jo3','j o3->o1d',    &
     411                             's-1',3,jo3_3d(1,1))
    495412           if (callstats) then
    496             call wstats(ngridmx,'jo3','j o3->o1d',
    497      $                       's-1',3,jo3_3d(1,1))
     413              call wstats(ngridmx,'jo3','j o3->o1d',       &
     414                          's-1',3,jo3_3d(1,1))
    498415           endif
    499416         end if ! of if (ngridmx.gt.1)
    500       endif ! of if (output)
    501 c
     417      end if ! of if (output)
     418
     419      return
    502420      end

  • trunk/LMDZ.MARS/libf/aeronomars/concentrations.F

    r370 r618  
    33      implicit none
    44
    5 c=======================================================================
    6 c CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
    7 c
    8 c mmean(ngridmx,nlayermx)       amu
    9 c cpnew(ngridmx,nlayermx)       J/kg/K
    10 c rnew(ngridmx,nlayermx)        J/kg/K
    11 c akknew(ngridmx,nlayermx)      coefficient of thermal concduction
    12 c
    13 c version: March 2011 - Franck Lefevre
    14 c=======================================================================
     5!=======================================================================
     6! CALCULATION OF MEAN MOLECULAR MASS, Cp, Akk and R
     7!
     8! mmean(ngridmx,nlayermx)       amu
     9! cpnew(ngridmx,nlayermx)       J/kg/K
     10! rnew(ngridmx,nlayermx)        J/kg/K
     11! akknew(ngridmx,nlayermx)      coefficient of thermal concduction
     12!
     13! version: April 2012 - Franck Lefevre
     14!=======================================================================
    1515
    16 c    Declarations
    17 c    ------------
     16!     declarations
    1817 
    1918#include "dimensions.h"
    … …  
    2625#include "conc.h"
    2726
    28 c    Input/Output
    29 c    ------------
     27!     input/output
    3028
    3129      real pplay(ngridmx,nlayermx)
    … …  
    3634      real ptimestep
    3735
    38 c    Local variables
    39 c    ---------------
     36!     local variables
    4037
    41       integer       :: l, ig, n, k
    42       integer, save :: gind(ncomp)
     38      integer       :: i, l, ig, iq
     39      integer, save :: nbq, niq(nqmx)
    4340      real          :: ni(nqmx), ntot
    44       real          :: zq(ngridmx,nlayermx,ncomp)
    45       real          :: zt(ngridmx,nlayermx)
    46       real, save    :: aki(ncomp)
    47       real, save    :: cpi(ncomp)
     41      real          :: zq(ngridmx, nlayermx, nqmx)
     42      real          :: zt(ngridmx, nlayermx)
     43      real, save    :: aki(nqmx)
     44      real, save    :: cpi(nqmx)
    4845
    4946      logical, save :: firstcall = .true.
    5047
    51 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    52 c     tracer numbering for the thermal conduction and
    53 c     specific heat coefficients
    54 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    55 
    56       integer,parameter :: i_co2  = 1
    57       integer,parameter :: i_co   = 2
    58       integer,parameter :: i_o    = 3
    59       integer,parameter :: i_o1d  = 4
    60       integer,parameter :: i_o2   = 5
    61       integer,parameter :: i_o3   = 6
    62       integer,parameter :: i_h    = 7
    63       integer,parameter :: i_h2   = 8
    64       integer,parameter :: i_oh   = 9
    65       integer,parameter :: i_ho2  = 10
    66       integer,parameter :: i_h2o2 = 11
    67       integer,parameter :: i_ch4  = 12
    68       integer,parameter :: i_n2   = 13
    69       integer,parameter :: i_ar   = 14
    70       integer,parameter :: i_h2o  = 15
    71 
    7248      if (firstcall) then
    7349
    74 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    75 c        initializations at first call:
    76 c        fill local array of tracer indexes
    77 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     50!        find index of chemical tracers to use
     51!        initialize thermal conductivity and specific heat coefficients
     52!        !? values are estimated
    7853
    79          gind(i_co2)   =  igcm_co2     ! co2
    80          gind(i_co)    =  igcm_co      ! co
    81          gind(i_o)     =  igcm_o       ! o
    82          gind(i_o1d)   =  igcm_o1d     ! o1d
    83          gind(i_o2)    =  igcm_o2      ! o2
    84          gind(i_o3)    =  igcm_o3      ! o3
    85          gind(i_h)     =  igcm_h       ! h
    86          gind(i_h2)    =  igcm_h2      ! h2
    87          gind(i_oh)    =  igcm_oh      ! oh
    88          gind(i_ho2)   =  igcm_ho2     ! ho2
    89          gind(i_h2o2)  =  igcm_h2o2    ! h2o2
    90          gind(i_ch4)   =  igcm_ch4     ! ch4
    91          gind(i_n2)    =  igcm_n2      ! n2
    92          gind(i_ar)    =  igcm_ar      ! ar
    93          gind(i_h2o)   =  igcm_h2o_vap ! h2o
     54         nbq = 0 ! to count number of tracers used in this subroutine
    9455
    95 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    96 c    Thermal conductivity and specific heat coefficients
    97 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     56         if (igcm_co2 /= 0) then
     57            nbq = nbq + 1
     58            niq(nbq) = igcm_co2
     59            aki(nbq) = 3.072e-4
     60            cpi(nbq) = 0.834e3
     61         end if
     62         if (igcm_co /= 0) then
     63            nbq = nbq + 1
     64            niq(nbq) = igcm_co
     65            aki(nbq) = 4.87e-4
     66            cpi(nbq) = 1.034e3
     67         end if
     68         if (igcm_o /= 0) then
     69            nbq = nbq + 1
     70            niq(nbq) = igcm_o
     71            aki(nbq) = 7.59e-4
     72            cpi(nbq) = 1.3e3
     73         end if
     74         if (igcm_o1d /= 0) then
     75            nbq = nbq + 1
     76            niq(nbq) = igcm_o1d
     77            aki(nbq) = 7.59e-4  !?
     78            cpi(nbq) = 1.3e3    !?
     79         end if
     80         if (igcm_o2 /= 0) then
     81            nbq = nbq + 1
     82            niq(nbq) = igcm_o2
     83            aki(nbq) = 5.68e-4
     84            cpi(nbq) = 0.9194e3
     85         end if
     86         if (igcm_o3 /= 0) then
     87            nbq = nbq + 1
     88            niq(nbq) = igcm_o3
     89            aki(nbq) = 3.00e-4  !?
     90            cpi(nbq) = 0.800e3  !?
     91         end if
     92         if (igcm_h /= 0) then
     93            nbq = nbq + 1
     94            niq(nbq) = igcm_h
     95            aki(nbq) = 0.0
     96            cpi(nbq) = 20.780e3
     97         end if
     98         if (igcm_h2 /= 0) then
     99            nbq = nbq + 1
     100            niq(nbq) = igcm_h2
     101            aki(nbq) = 36.314e-4
     102            cpi(nbq) = 14.266e3
     103         end if
     104         if (igcm_oh /= 0) then
     105            nbq = nbq + 1
     106            niq(nbq) = igcm_oh
     107            aki(nbq)  = 7.00e-4 !?
     108            cpi(nbq)  = 1.045e3
     109         end if
     110         if (igcm_ho2 /= 0) then
     111            nbq = nbq + 1
     112            niq(nbq) = igcm_ho2
     113            aki(nbq) = 0.0
     114            cpi(nbq) = 1.065e3  !?
     115         end if
     116         if (igcm_n2 /= 0) then
     117            nbq = nbq + 1
     118            niq(nbq) = igcm_n2
     119            aki(nbq) = 5.6e-4
     120            cpi(nbq) = 1.034e3
     121         end if
     122         if (igcm_ar /= 0) then
     123            nbq = nbq + 1
     124            niq(nbq) = igcm_ar
     125            aki(nbq) = 0.0      !?
     126            cpi(nbq) = 1.000e3  !?
     127         end if
     128         if (igcm_h2o_vap /= 0) then
     129            nbq = nbq + 1
     130            niq(nbq) = igcm_h2o_vap
     131            aki(nbq) = 0.0
     132            cpi(nbq) = 1.870e3
     133         end if
    98134
    99          aki(i_co2)   = 3.072e-4
    100          aki(i_co)    = 4.87e-4
    101          aki(i_o)     = 7.59e-4
    102          aki(i_o1d)   = 7.59e-4    !?
    103          aki(i_o2)    = 5.68e-4
    104          aki(i_o3)    = 3.00e-4    !?
    105          aki(i_h)     = 0.0
    106          aki(i_h2)    = 36.314e-4
    107          aki(i_oh)    = 7.00e-4    !?
    108          aki(i_ho2)   = 0.0
    109          aki(i_h2o2)  = 0.0
    110          aki(i_ch4)   = 0.0        !?
    111          aki(i_n2)    = 5.6e-4
    112          aki(i_ar)    = 0.0        !?
    113          aki(i_h2o)   = 0.0
     135         firstcall = .false.
    114136
    115          cpi(i_co2)   = 0.834e3
    116          cpi(i_co)    = 1.034e3
    117          cpi(i_o)     = 1.3e3
    118          cpi(i_o1d)   = 1.3e3    !?
    119          cpi(i_o2)    = 0.9194e3
    120          cpi(i_o3)    = 0.800e3  !?
    121          cpi(i_h)     = 20.780e3
    122          cpi(i_h2)    = 14.266e3
    123          cpi(i_oh)    = 1.045e3
    124          cpi(i_ho2)   = 1.065e3  !?
    125          cpi(i_h2o2)  = 1.000e3  !?
    126          cpi(i_ch4)   = 1.000e3  !?
    127          cpi(i_n2)    = 1.034e3
    128          cpi(i_ar)    = 1.000e3  !?
    129          cpi(i_h2o)   = 1.870e3
     137      end if ! if (firstcall)
    130138
    131          firstcall=.false.
     139!     update temperature
    132140
    133       end if ! of if (firstcall)
    134 c
    135 c     initializations
    136 c
    137       mmean(:,:)  = 0.
    138       cpnew(:,:)  = 0.
    139       akknew(:,:) = 0.
    140 c
    141 c     update temperature
    142 c
    143141      do l = 1,nlayermx
    144142         do ig = 1,ngridmx
    … …  
    146144         end do
    147145      end do
    148 c
    149 c     update tracers
    150 c
     146
     147!     update tracers
     148
    151149      do l = 1,nlayermx
    152150         do ig = 1,ngridmx
    153             do n = 1,ncomp
    154                zq(ig,l,n) = max(1.e-30, pq(ig,l,gind(n))
    155      $                                + pdq(ig,l,gind(n))*ptimestep)
     151            do i = 1,nbq
     152               iq = niq(i)
     153               zq(ig,l,iq) = max(1.e-30, pq(ig,l,iq)
     154     $                                 + pdq(ig,l,iq)*ptimestep)
    156155            end do
    157156         end do
    158157      end do
    159 c
    160 c     mmean : mean molecular mass
    161 c     rnew  : specific gas constant
    162 c
     158
     159!     mmean : mean molecular mass
     160!     rnew  : specific gas constant
     161
     162      mmean(:,:)  = 0.
     163
    163164      do l = 1,nlayermx
    164165         do ig = 1,ngridmx
    165             do n = 1, ncomp
    166                mmean(ig,l) = mmean(ig,l) + zq(ig,l,n)/mmol(gind(n))
     166            do i = 1,nbq
     167               iq = niq(i)
     168               mmean(ig,l) = mmean(ig,l) + zq(ig,l,iq)/mmol(iq)
    167169            end do
    168170            mmean(ig,l) = 1./mmean(ig,l)
    … …  
    170172         end do
    171173      end do
    172 c
    173 c     cpnew  : specicic heat
    174 c     akknew : thermal conductivity cofficient
    175 c     
     174
     175!     cpnew  : specicic heat
     176!     akknew : thermal conductivity cofficient
     177     
     178      cpnew(:,:)  = 0.
     179      akknew(:,:) = 0.
     180
    176181      do l = 1,nlayermx
    177182         do ig = 1,ngridmx
    178183            ntot = pplay(ig,l)/(1.381e-23*zt(ig,l))*1.e-6  ! in #/cm3
    179             do n = 1,ncomp
    180                ni(n) = ntot*zq(ig,l,n)*mmean(ig,l)/mmol(gind(n))
    181                cpnew(ig,l) = cpnew(ig,l) + ni(n)*cpi(n)
    182                akknew(ig,l) = akknew(ig,l) + ni(n)*aki(n)
     184            do i = 1,nbq
     185               iq = niq(i)
     186               ni(iq) = ntot*zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
     187               cpnew(ig,l) = cpnew(ig,l) + ni(iq)*cpi(iq)
     188               akknew(ig,l) = akknew(ig,l) + ni(iq)*aki(iq)
    183189            end do
    184190            cpnew(ig,l) = cpnew(ig,l)/ntot
    185191            akknew(ig,l) = akknew(ig,l)/ntot
    186192         end do
    187 c        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
     193!        print*, l, mmean(1,l), cpnew(1,l), rnew(1,l)
    188194      end do
    189195

  • trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90

    r616 r618  
    1       subroutine inichim_newstart(pq,ps,sig,nq,latfi,lonfi,airefi,
    2      $ thermo,qsurf)
     1      subroutine inichim_newstart(pq, qsurf, ps, flagh2o, flagthermo)
    32
    43      implicit none
    54
    6 c=======================================================================
    7 c
    8 c   subject:
    9 c   --------
    10 c
    11 c  Initialization of the composition for use in the building of a new start file
    12 c  used by program newstart.F
    13 c  also used by program testphys1d.F
    14 c
    15 c   Author:   Sebastien Lebonnois (08/11/2002)
    16 c   -------
    17 c   Revised 07/2003 by Monica Angelats-i-Coll to use input files
    18 c   Modified 10/2008 Identify tracers by their names Ehouarn Millour
    19 c   Modified 11/2011 Addition of methane (Franck Lefevre)
    20 c
    21 c   Arguments:
    22 c   ----------
    23 c
    24 c    pq(iip1,jjp1,llm,nqmx)  Advected fields, ie chemical species here
    25 c    sig = sigma                vertical coordinate (interface of layers)
    26 c    nq: number of tracers to initialize (used to evaluate if only
    27 c        chemistry species are to be initialized or if water vapour
    28 c        and water ice should also be initialized.
    29 c        NB: number of chemistry tracers is ncomp (set in chimiedata.h)
    30 c=======================================================================
    31 
    32 c    Declarations :
    33 c    --------------
     5!=======================================================================
     6!
     7!  Purpose:
     8!  --------
     9!
     10!  Initialization of the chemistry for use in the building of a new start file
     11!  used by program newstart.F
     12!  also used by program testphys1d.F
     13!
     14!  Authors:
     15!  -------
     16!  Initial version 11/2002 by Sebastien Lebonnois
     17!  Revised 07/2003 by Monica Angelats-i-Coll to use input files
     18!  Modified 10/2008 Identify tracers by their names Ehouarn Millour
     19!  Modified 11/2011 Addition of methane Franck Lefevre
     20!  Rewritten 04/2012 Franck Lefevre
     21!
     22!  Arguments:
     23!  ----------
     24!
     25!  pq(iip1,jjp1,llm,nqmx)  Advected fields, ie chemical species here
     26!  qsurf(ngridmx,nqmx)     Amount of tracer on the surface (kg/m2)
     27!  ps(iip1,jjp1)           Surface pressure (Pa)
     28!  flagh2o                 flag for initialisation of h2o (1: yes / 0: no)
     29!  flagthermo              flag for initialisation of thermosphere only (1: yes / 0: no)
     30!
     31!=======================================================================
    3432
    3533#include "dimensions.h"
    3634#include "dimphys.h"
    3735#include "paramet.h"
    38 #include "chimiedata.h"
    3936#include "tracer.h"
    40 #include "comcstfi.h"
    41 #include "comdiurn.h"
     37#include "comvert.h"
    4238#include "callkeys.h"
    43 #include "temps.h"
    4439#include "datafile.h"
    4540
    46 c Arguments :
    47 c -----------
    48 
    49       real    ps(iip1,jjp1)
    50       real    pq(iip1,jjp1,llm,nqmx)     ! Advected fields, ie chemical species
    51       real    sig(llm+1)                 ! vertical coordinate (interface of layers)
    52       integer nq            !  =nqmx
    53                             ! or =nqmx-1 if H2O kept
    54                             ! or =nqmx-2 if H2O and ice kept
    55       REAL latfi(ngridmx),lonfi(ngridmx),airefi(ngridmx)
    56       integer thermo ! flag for thermosphere initialisation only
    57       real :: qsurf(ngridmx,nqmx) ! surface tracers
    58 
    59 c Local variables :
    60 c -----------------
    61 
    62       integer  iq,i,j,l,n, nspecini,inil,iqmax,iiq
    63       INTEGER aa(nqmx)
    64       REAL qtot(iip1,jjp1,llm)
    65       real densconc,zgcm,zfile(252),vmrint,nt, mmean
    66       real nxf(252),zfilemin(252),sigfile(252)
    67       real vmrf(252,14),nf(252,14)
    68       real tfile(252),zzfile(252)
    69       real ni(14),niold(14)
    70       real mmolf(14)    ! molecular mass in amu (species in files)
    71       data mmolf/44.,40.,28.,32.,28.,16.,2.,   ! majors
    72      &           1.,17.,33.,18.,34.,16.,48./   ! minors
    73       character*3 tmp ! temporary variable
    74       integer ierr
    75 
    76       logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
    77       character(len=20) :: txt ! to store some text
    78       integer :: nqchem_start,count
    79       integer :: nqchem(nqmx) ! local chemistry tracer index array
    80       integer :: nbqchem ! total number of chemical species (water included)
    81 
    82 c-----------------------------------------------------------------------
    83 c    Input/Output
    84 c    ------------
    85 ! read in 'callphys.def'
    86        call inichim_readcallphys()
    87 
    88       ! check if tracers have 'old' names
    89       oldnames=.false.
    90       count=0
    91       do iq=1,nqmx
    92         txt=" "
    93         write(txt,'(a1,i2.2)') 'q',iq
    94         if (txt.eq.noms(iq)) then
    95           count=count+1
    96         endif
    97       enddo ! of do iq=1,nqmx
    98 
    99       if (count.eq.nqmx) then
    100         write(*,*) "inichim_newstart: tracers seem to follow old ",
    101      &             "naming convention (q01,q02,...)"
    102         write(*,*) "   => will work for now ... "
    103         write(*,*) "      but you should rename them"
    104         oldnames=.true.
    105       endif
    106 
    107       ! copy/set tracer names
    108       if (oldnames) then ! old names (derived from iq & option values)
    109         ! 1. dust:
    110         if (dustbin.ne.0) then ! transported dust
    111           do iq=1,dustbin
    112             txt=" "
    113             write(txt,'(a4,i2.2)') 'dust',iq
    114             noms(iq)=txt
    115             mmol(iq)=100.
    116           enddo
    117           ! special case if doubleq
    118           if (doubleq) then
    119             if (dustbin.ne.2) then
    120               write(*,*) 'initracer: error expected dustbin=2'
    121             else
    122 !              noms(1)='dust_mass'   ! dust mass mixing ratio
    123 !              noms(2)='dust_number' ! dust number mixing ratio
    124 ! same as above, but avoid explict possible out-of-bounds on array
    125                noms(1)='dust_mass'   ! dust mass mixing ratio
    126                do iq=2,2
    127                noms(iq)='dust_number' ! dust number mixing ratio
    128                enddo
    129             endif
    130           endif
    131         endif
    132         ! 2. water & ice
    133         if (water) then
    134           noms(nqmx)='h2o_vap'
    135           mmol(nqmx)=18.
    136 !            noms(nqmx-1)='h2o_ice'
    137 !            mmol(nqmx-1)=18.
    138           !"loop" to avoid potential out-of-bounds in array
    139           do iq=nqmx-1,nqmx-1
    140             noms(iq)='h2o_ice'
    141             mmol(iq)=18.
    142           enddo
    143         endif
    144         ! 3. Chemistry
    145         if (photochem .or. callthermos) then
    146           if (doubleq) then
    147             nqchem_start=3
    148           else
    149             nqchem_start=dustbin+1
    150           end if
    151         endif ! of if (photochem .or. callthermos)
    152         do iq=nqchem_start,nqchem_start+ncomp-2
    153           noms(iq)=nomchem(iq-nqchem_start+1)
    154           mmol(iq)=mmolchem(iq-nqchem_start+1)
    155         enddo
    156         ! 4. Other tracers ????
    157         if ((dustbin.eq.0).and.(.not.water)) then
    158           noms(1)='co2'
    159           mmol(1)=44
    160           if (nqmx.eq.2) then
    161             noms(nqmx)='Ar_N2'
    162             mmol(nqmx)=30
    163           endif
    164         endif
     41! inputs :
     42
     43      real,intent(in) :: ps(iip1,jjp1)            ! surface pressure in the gcm (Pa)   
     44      integer,intent(in) :: flagh2o               ! flag for h2o initialisation
     45      integer,intent(in) :: flagthermo            ! flag for thermosphere initialisation only
     46
     47! outputs :
     48
     49      real,intent(out) :: pq(iip1,jjp1,llm,nqmx)  ! advected fields, ie chemical species
     50      real,intent(out) :: qsurf(ngridmx,nqmx)     ! surface values (kg/m2) of tracers
     51
     52! local :
     53
     54      integer :: iq, i, j, l, n, nbqchem
     55      integer :: count, ierr, dummy
     56      real    :: mmean(iip1,jjp1,llm)             ! mean molecular mass (g)
     57      real    :: pgcm                             ! pressure at each layer in the gcm (Pa)
     58
     59      integer, parameter         :: nalt = 252    ! number of altitudes in the initialization files
     60      integer, parameter         :: nspe = 14     ! number of species in the initialization files
     61      integer, dimension(nspe)   :: niq           ! local index of species in initialization files
     62      real, dimension(nalt)      :: tinit, zzfile ! temperature in initialization files
     63      real, dimension(nalt)      :: pinit         ! pressure in initialization files
     64      real, dimension(nalt)      :: densinit      ! total number density in initialization files
     65      real, dimension(nalt,nspe) :: vmrinit       ! mixing ratios in initialization files
     66      real, dimension(nspe)      :: vmrint        ! mixing ratio interpolated onto the gcm vertical grid
     67      real                       :: vmr
     68
     69      character(len=20)          :: txt           ! to store some text
     70
     71! 1. identify tracers by their names: (and set corresponding values of mmol)
     72
     73! 1.1 initialize tracer indexes to zero:
     74
     75      do iq = 1,nqmx
     76        igcm_dustbin(iq) = 0
     77      end do
     78
     79      igcm_dust_mass   = 0
     80      igcm_dust_number = 0
     81      igcm_ccn_mass    = 0
     82      igcm_ccn_number  = 0
     83      igcm_h2o_vap     = 0
     84      igcm_h2o_ice     = 0
     85      igcm_co2         = 0
     86      igcm_co          = 0
     87      igcm_o           = 0
     88      igcm_o1d         = 0
     89      igcm_o2          = 0
     90      igcm_o3          = 0
     91      igcm_h           = 0
     92      igcm_h2          = 0
     93      igcm_oh          = 0
     94      igcm_ho2         = 0
     95      igcm_h2o2        = 0
     96      igcm_ch4         = 0
     97      igcm_n2          = 0
     98      igcm_ar          = 0
     99      igcm_ar_n2       = 0
     100      igcm_co2plus     = 0
     101      igcm_oplus       = 0
     102      igcm_o2plus      = 0
     103      igcm_coplus      = 0
     104      igcm_cplus       = 0
     105      igcm_nplus       = 0
     106      igcm_noplus      = 0
     107      igcm_n2plus      = 0
     108      igcm_hplus       = 0
     109      igcm_elec        = 0
     110
     111! 1.2 find dust tracers
     112
     113      count = 0
     114
     115      if (dustbin > 0) then
     116         do iq = 1,nqmx
     117            txt = " "
     118            write(txt,'(a4,i2.2)') 'dust', count + 1
     119            if (noms(iq) == txt) then
     120               count = count + 1
     121               igcm_dustbin(count) = iq
     122               mmol(iq) = 100.
     123            end if
     124         end do !do iq=1,nqmx
     125      end if ! of if (dustbin.gt.0)
     126
     127      if (doubleq) then
     128         do iq = 1,nqmx
     129            if (noms(iq) == "dust_mass") then
     130               igcm_dust_mass = iq
     131               count = count + 1
     132            end if
     133            if (noms(iq) == "dust_number") then
     134               igcm_dust_number = iq
     135               count = count + 1
     136            end if
     137         end do
     138      end if ! of if (doubleq)
     139
     140      if (scavenging) then
     141         do iq = 1,nqmx
     142            if (noms(iq) == "ccn_mass") then
     143               igcm_ccn_mass = iq
     144               count = count + 1
     145            end if
     146            if (noms(iq) == "ccn_number") then
     147               igcm_ccn_number = iq
     148               count = count + 1
     149            end if
     150         end do
     151      end if ! of if (scavenging)
     152
     153      if (submicron) then
     154         do iq=1,nqmx
     155            if (noms(iq) == "dust_submicron") then
     156               igcm_dust_submicron = iq
     157               mmol(iq) = 100.
     158               count = count + 1
     159            end if
     160         end do
     161      end if ! of if (submicron)
     162
     163! 1.3 find chemistry and water tracers
     164
     165      nbqchem = 0
     166      do iq = 1,nqmx
     167         if (noms(iq) == "co2") then
     168            igcm_co2 = iq
     169            mmol(igcm_co2) = 44.
     170            count = count + 1
     171            nbqchem = nbqchem + 1
     172        end if
     173        if (noms(iq) == "co") then
     174           igcm_co = iq
     175           mmol(igcm_co) = 28.
     176           count = count + 1
     177           nbqchem = nbqchem + 1
     178        end if
     179        if (noms(iq) == "o") then
     180           igcm_o = iq
     181           mmol(igcm_o) = 16.
     182           count = count + 1
     183           nbqchem = nbqchem + 1
     184        end if
     185        if (noms(iq) == "o1d") then
     186           igcm_o1d = iq
     187           mmol(igcm_o1d) = 16.
     188           count = count + 1
     189           nbqchem = nbqchem + 1
     190        end if
     191        if (noms(iq) == "o2") then
     192           igcm_o2 = iq
     193           mmol(igcm_o2) = 32.
     194           count = count + 1
     195           nbqchem = nbqchem + 1
     196        end if
     197        if (noms(iq) == "o3") then
     198           igcm_o3 = iq
     199           mmol(igcm_o3) = 48.
     200           count = count + 1
     201           nbqchem = nbqchem + 1
     202        end if
     203        if (noms(iq) == "h") then
     204           igcm_h = iq
     205           mmol(igcm_h) = 1.
     206           count = count + 1
     207           nbqchem = nbqchem + 1
     208        end if
     209        if (noms(iq) == "h2") then
     210           igcm_h2 = iq
     211           mmol(igcm_h2) = 2.
     212           count = count + 1
     213           nbqchem = nbqchem + 1
     214        end if
     215        if (noms(iq) == "oh") then
     216           igcm_oh = iq
     217           mmol(igcm_oh) = 17.
     218           count = count + 1
     219           nbqchem = nbqchem + 1
     220        end if
     221        if (noms(iq) == "ho2") then
     222           igcm_ho2 = iq
     223           mmol(igcm_ho2) = 33.
     224           count = count + 1
     225           nbqchem = nbqchem + 1
     226        end if
     227        if (noms(iq) == "h2o2") then
     228           igcm_h2o2 = iq
     229           mmol(igcm_h2o2) = 34.
     230           count = count + 1
     231           nbqchem = nbqchem + 1
     232        end if
     233        if (noms(iq) == "ch4") then
     234           igcm_ch4 = iq
     235           mmol(igcm_ch4) = 16.
     236           count = count + 1
     237           nbqchem = nbqchem + 1
     238        end if
     239        if (noms(iq) == "n2") then
     240           igcm_n2 = iq
     241           mmol(igcm_n2) = 28.
     242           count = count + 1
     243           nbqchem = nbqchem + 1
     244        end if
     245        if (noms(iq) == "n") then
     246           igcm_n = iq
     247           mmol(igcm_n) = 14.
     248           count = count + 1
     249           nbqchem = nbqchem + 1
     250        end if
     251        if (noms(iq) == "n2d") then
     252           igcm_n2d = iq
     253           mmol(igcm_n2d) = 14.
     254           count = count + 1
     255           nbqchem = nbqchem + 1
     256        end if
     257        if (noms(iq) == "no") then
     258           igcm_no = iq
     259           mmol(igcm_no) = 30.
     260           count = count + 1
     261           nbqchem = nbqchem + 1
     262        end if
     263        if (noms(iq) == "no2") then
     264           igcm_no2 = iq
     265           mmol(igcm_no2) = 46.
     266           count = count + 1
     267           nbqchem = nbqchem + 1
     268        end if
     269        if (noms(iq) == "ar") then
     270           igcm_ar = iq
     271           mmol(igcm_ar) = 40.
     272           count = count + 1
     273           nbqchem = nbqchem + 1
     274        end if
     275        if (noms(iq) == "h2o_vap") then
     276           igcm_h2o_vap = iq
     277           mmol(igcm_h2o_vap) = 18.
     278           count = count + 1
     279           nbqchem = nbqchem + 1
     280        end if
     281        if (noms(iq) == "h2o_ice") then
     282           igcm_h2o_ice = iq
     283           mmol(igcm_h2o_ice) = 18.
     284           count = count + 1
     285           nbqchem = nbqchem + 1
     286        end if
     287
     288! 1.4 find ions
     289
     290        if (noms(iq) == "co2plus") then
     291           igcm_co2plus = iq
     292           mmol(igcm_co2plus) = 44.
     293           count = count + 1
     294           nbqchem = nbqchem + 1
     295        end if
     296        if (noms(iq) == "oplus") then
     297           igcm_oplus = iq
     298           mmol(igcm_oplus) = 16.
     299           count = count + 1
     300           nbqchem = nbqchem + 1
     301        end if
     302        if (noms(iq) == "o2plus") then
     303           igcm_o2plus = iq
     304           mmol(igcm_o2plus) = 32.
     305           count = count + 1
     306           nbqchem = nbqchem + 1
     307        end if
     308        if (noms(iq) == "coplus") then
     309           igcm_coplus = iq
     310           mmol(igcm_coplus) = 28.
     311           count = count + 1
     312           nbqchem = nbqchem + 1
     313        end if
     314        if (noms(iq) == "cplus") then
     315           igcm_cplus = iq
     316           mmol(igcm_cplus) = 12.
     317           count = count + 1
     318           nbqchem = nbqchem + 1
     319        end if
     320        if (noms(iq) == "nplus") then
     321           igcm_nplus = iq
     322           mmol(igcm_nplus) = 14.
     323           count = count + 1
     324           nbqchem = nbqchem + 1
     325        end if
     326        if (noms(iq) == "noplus") then
     327           igcm_noplus = iq
     328           mmol(igcm_noplus) = 30.
     329           count = count + 1
     330           nbqchem = nbqchem + 1
     331        end if
     332        if (noms(iq) == "n2plus") then
     333           igcm_n2plus = iq
     334           mmol(igcm_n2plus) = 28.
     335           count = count + 1
     336           nbqchem = nbqchem + 1
     337        end if
     338        if (noms(iq) == "hplus") then
     339           igcm_hplus = iq
     340           mmol(igcm_hplus) = 1.
     341           count = count + 1
     342           nbqchem = nbqchem + 1
     343        end if
     344        if (noms(iq) == "elec") then
     345           igcm_elec = iq
     346           mmol(igcm_elec) = 1./1822.89
     347           count = count + 1
     348        end if
     349
     350! 1.5 find idealized non-condensible tracer
     351
     352        if (noms(iq) == "Ar_N2") then
     353           igcm_ar_n2 = iq
     354           mmol(igcm_ar_n2) = 30.
     355           count = count + 1
     356        end if
     357
     358      end do ! of do iq=1,nqmx
     359     
     360! 1.6 check that we identified all tracers:
     361
     362      if (count /= nqmx) then
     363         write(*,*) "inichim_newstart: found only ",count," tracers"
     364         write(*,*) "                  expected ",nqmx
     365         do iq = 1,count
     366            write(*,*) '      ', iq, ' ', trim(noms(iq))
     367         end do
     368         stop
    165369      else
    166         ! noms(:) already contain tracer names
    167         ! mmol(:) array is initialized later (see below)
    168       endif ! of if (oldnames)
    169 
    170       ! special modification when using 'old' tracers:
    171       ! qsurf(nqmx) was h2o ice whereas q(nqmx) was water vapour
    172       ! and (if iceparty) q(nqmx-1) was null whereas q(nqmx-1) was water ice
    173       if (oldnames.and.water) then
    174         write(*,*)'inichim_newstart: moving surface water ice to index '
    175      &            ,nqmx-1
    176         ! "loop" to avoid potential out-of-bounds on arrays
    177         do iq=nqmx-1,nqmx-1
    178         qsurf(1:ngridmx,iq)=qsurf(1:ngridmx,iq+1)
    179         enddo
    180         qsurf(1:ngridmx,nqmx)=0
    181       endif
    182 
    183 c Dimension test:
    184 c ---------------
    185 
    186        if (water) then
    187 !         if ((nqchem_min+ncomp+1).ne.nqmx) then
    188          if ((dustbin+ncomp+2).ne.nqmx) then
    189 !            print*,'********* Dimension problem! ********'
    190 !            print*,"nqchem_min+ncomp+1).ne.nqmx"
    191             print*,'inichim_newstart: tracer dimension problem:'
    192             print*,"(dustbin+ncomp+2).ne.nqmx"
    193             print*,"ncomp: ",ncomp
    194 !            print*,"nqchem_min: ",nqchem_min
    195             print*,"nqmx: ",nqmx
    196             print*,'Change ncomp in chimiedata.h'
    197             STOP
    198          endif
    199        else
    200 !         if ((nqchem_min+ncomp).ne.nqmx) then
    201           if ((dustbin+ncomp+1).ne.nqmx) then
    202 !            print*,'********* Dimension problem! ********'
    203 !            print*,"nqchem_min+ncomp).ne.nqmx"
    204             print*,'initracer: tracer dimension problem:'
    205             print*,"(dustbin+ncomp+1).ne.nqmx"
    206             print*,"ncomp: ",ncomp
    207 !            print*,"nqchem_min: ",nqchem_min
    208             print*,"nqmx: ",nqmx
    209             print*,'Change ncomp in chimiedata.h'
    210             STOP
    211          endif
    212        endif
    213 
    214 c noms and mmol vectors:
    215 c ----------------------
    216 !
    217 !       if (iceparty) then
    218 !          do iq=nqchem_min,nqmx-2
    219 !            noms(iq) = nomchem(iq-nqchem_min+1)
    220 !            mmol(iq) = mmolchem(iq-nqchem_min+1)
    221 !          enddo
    222 !          noms(nqmx-1) = 'ice'
    223 !          mmol(nqmx-1) = 18.
    224 !          noms(nqmx)   = 'h2o'
    225 !          mmol(nqmx)   = 18.
    226 !       else
    227 !          do iq=nqchem_min,nqmx-1
    228 !            noms(iq) = nomchem(iq-nqchem_min+1)
    229 !            mmol(iq) = mmolchem(iq-nqchem_min+1)
    230 !          enddo
    231 !          noms(nqmx)   = 'h2o'
    232 !          mmol(nqmx)   = 18.
    233 !       endif
    234 !       if (nqchem_min.gt.1) then
    235 !          do iq=1,nqchem_min-1
    236 !            noms(iq) = 'dust'
    237 !            mmol(iq) = 100.
    238 !          enddo
    239 !       endif
    240 
    241 
    242 ! Identify tracers by their names: (and set corresponding values of mmol)
    243       ! 0. initialize tracer indexes to zero:
    244       do iq=1,nqmx
    245         igcm_dustbin(iq)=0
    246       enddo
    247       igcm_dust_mass=0
    248       igcm_dust_number=0
    249       igcm_h2o_vap=0
    250       igcm_h2o_ice=0
    251       igcm_co2=0
    252       igcm_co=0
    253       igcm_o=0
    254       igcm_o1d=0
    255       igcm_o2=0
    256       igcm_o3=0
    257       igcm_h=0
    258       igcm_h2=0
    259       igcm_oh=0
    260       igcm_ho2=0
    261       igcm_h2o2=0
    262       igcm_ch4=0
    263       igcm_n2=0
    264       igcm_ar=0
    265       igcm_ar_n2=0
    266 
    267       ! 1. find dust tracers
    268       count=0
    269       if (dustbin.gt.0) then
    270         do iq=1,nqmx
    271           txt=" "
    272           write(txt,'(a4,i2.2)')'dust',count+1
    273           if (noms(iq).eq.txt) then
    274             count=count+1
    275             igcm_dustbin(count)=iq
    276             mmol(iq)=100.
    277           endif
    278         enddo !do iq=1,nqmx
    279       endif ! of if (dustbin.gt.0)
    280       if (doubleq) then
    281         do iq=1,nqmx
    282           if (noms(iq).eq."dust_mass") then
    283             igcm_dust_mass=iq
    284             count=count+1
    285           endif
    286           if (noms(iq).eq."dust_number") then
    287             igcm_dust_number=iq
    288             count=count+1
    289           endif
    290         enddo
    291       endif ! of if (doubleq)
    292       ! 2. find chemistry and water tracers
    293       nbqchem=0
    294       do iq=1,nqmx
    295         if (noms(iq).eq."co2") then
    296           igcm_co2=iq
    297           mmol(igcm_co2)=44.
    298           count=count+1
    299           nbqchem=nbqchem+1
    300         endif
    301         if (noms(iq).eq."co") then
    302           igcm_co=iq
    303           mmol(igcm_co)=28.
    304           count=count+1
    305           nbqchem=nbqchem+1
    306         endif
    307         if (noms(iq).eq."o") then
    308           igcm_o=iq
    309           mmol(igcm_o)=16.
    310           count=count+1
    311           nbqchem=nbqchem+1
    312         endif
    313         if (noms(iq).eq."o1d") then
    314           igcm_o1d=iq
    315           mmol(igcm_o1d)=16.
    316           count=count+1
    317           nbqchem=nbqchem+1
    318         endif
    319         if (noms(iq).eq."o2") then
    320           igcm_o2=iq
    321           mmol(igcm_o2)=32.
    322           count=count+1
    323           nbqchem=nbqchem+1
    324         endif
    325         if (noms(iq).eq."o3") then
    326           igcm_o3=iq
    327           mmol(igcm_o3)=48.
    328           count=count+1
    329           nbqchem=nbqchem+1
    330         endif
    331         if (noms(iq).eq."h") then
    332           igcm_h=iq
    333           mmol(igcm_h)=1.
    334           count=count+1
    335           nbqchem=nbqchem+1
    336         endif
    337         if (noms(iq).eq."h2") then
    338           igcm_h2=iq
    339           mmol(igcm_h2)=2.
    340           count=count+1
    341           nbqchem=nbqchem+1
    342         endif
    343         if (noms(iq).eq."oh") then
    344           igcm_oh=iq
    345           mmol(igcm_oh)=17.
    346           count=count+1
    347           nbqchem=nbqchem+1
    348         endif
    349         if (noms(iq).eq."ho2") then
    350           igcm_ho2=iq
    351           mmol(igcm_ho2)=33.
    352           count=count+1
    353           nbqchem=nbqchem+1
    354         endif
    355         if (noms(iq).eq."h2o2") then
    356           igcm_h2o2=iq
    357           mmol(igcm_h2o2)=34.
    358           count=count+1
    359           nbqchem=nbqchem+1
    360         endif
    361         if (noms(iq).eq."ch4") then
    362           igcm_ch4=iq
    363           mmol(igcm_ch4)=16.
    364           count=count+1
    365           nbqchem=nbqchem+1
    366         endif
    367         if (noms(iq).eq."n2") then
    368           igcm_n2=iq
    369           mmol(igcm_n2)=28.
    370           count=count+1
    371           nbqchem=nbqchem+1
    372         endif
    373         if (noms(iq).eq."ar") then
    374           igcm_ar=iq
    375           mmol(igcm_ar)=40.
    376           count=count+1
    377           nbqchem=nbqchem+1
    378         endif
    379         if (noms(iq).eq."h2o_vap") then
    380           igcm_h2o_vap=iq
    381           mmol(igcm_h2o_vap)=18.
    382           count=count+1
    383           nbqchem=nbqchem+1
    384         endif
    385         if (noms(iq).eq."h2o_ice") then
    386           igcm_h2o_ice=iq
    387           mmol(igcm_h2o_ice)=18.
    388           count=count+1
    389           nbqchem=nbqchem+1
    390         endif
    391         ! Other stuff: e.g. for simulations using co2 + neutral gaz
    392         if (noms(iq).eq."Ar_N2") then
    393           igcm_ar_n2=iq
    394           mmol(igcm_ar_n2)=30.
    395           count=count+1
    396         endif
    397       enddo ! of do iq=1,nqmx
    398      
    399       ! check that we identified all tracers:
    400       if (count.ne.nqmx) then
    401         write(*,*) "inichim_newstart: found only ",count," tracers"
    402         write(*,*) "               expected ",nqmx
    403         do iq=1,count
    404           write(*,*)'      ',iq,' ',trim(noms(iq))
    405         enddo
    406         stop
    407       else
    408        write(*,*)"inichim_newstart: found all expected tracers, namely:"
    409         do iq=1,nqmx
    410           write(*,*)'      ',iq,' ',trim(noms(iq))
     370         write(*,*) "inichim_newstart: found all expected tracers"
     371         do iq = 1,nqmx
     372            write(*,*) '      ', iq, ' ', trim(noms(iq))
     373         end do
     374      end if
     375
     376! 2. load in chemistry data for initialization:
     377
     378! order of major species in initialization file:
     379!
     380!    1: co2
     381!    2: ar
     382!    3: n2 
     383!    4: o2 
     384!    5: co 
     385!    6: o   
     386!    7: h2
     387!
     388! order of minor species in initialization file:
     389!
     390!    1: h 
     391!    2: oh
     392!    3: ho2
     393!    4: h2o
     394!    5: h2o2
     395!    6: o1d
     396!    7: o3
     397
     398! major species:
     399
     400      niq(1) = igcm_co2
     401      niq(2) = igcm_ar
     402      niq(3) = igcm_n2
     403      niq(4) = igcm_o2
     404      niq(5) = igcm_co
     405      niq(6) = igcm_o
     406      niq(7) = igcm_h2
     407
     408! minor species:
     409
     410      niq(8)  = igcm_h
     411      niq(9)  = igcm_oh
     412      niq(10) = igcm_ho2
     413      niq(11) = igcm_h2o_vap
     414      niq(12) = igcm_h2o2
     415      niq(13) = igcm_o1d
     416      niq(14) = igcm_o3
     417
     418! 2.1 open initialization files
     419
     420      open(210, iostat=ierr,file=trim(datafile)// &
     421                            '/atmosfera_LMD_may.dat')
     422      if (ierr /= 0) then
     423         write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
     424         write(*,*)'(in aeronomars/inichim.F)'
     425         write(*,*)'It should be in :', trim(datafile),'/'
     426         write(*,*)'1) You can change this path in callphys.def with'
     427         write(*,*)'   datadir=/path/to/datafiles/'
     428         write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
     429         write(*,*)'   can be obtained online on:'
     430         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
     431         stop
     432      end if
     433      open(220, iostat=ierr,file=trim(datafile)// &
     434                            '/atmosfera_LMD_min.dat')
     435      if (ierr /= 0) then
     436         write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
     437         write(*,*)'(in aeronomars/inichim.F)'
     438         write(*,*)'It should be in :', trim(datafile),'/'
     439         write(*,*)'1) You can change this path in callphys.def with'
     440         write(*,*)'   datadir=/path/to/datafiles/'
     441         write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
     442         write(*,*)'   can be obtained online on:'
     443         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
     444         stop
     445      end if
     446
     447! 2.2 read initialization files
     448
     449! major species
     450
     451      read(210,*)
     452      do l = 1,nalt
     453         read(210,*) dummy, tinit(l), pinit(l), densinit(l), &
     454                     (vmrinit(l,n), n = 1,7)
     455         pinit(l) = pinit(l)*100.              ! conversion in Pa
     456         pinit(l) = log(pinit(l))              ! for the vertical interpolation
     457      end do
     458      close(210)
     459
     460! minor species
     461
     462      read(220,*)
     463      do l = 1,nalt
     464         read(220,*) dummy, (vmrinit(l,n), n = 8,14)
     465      end do
     466      close(220)
     467
     468! 3. initialization of tracers
     469
     470      do i = 1,iip1
     471         do j = 1,jjp1
     472            do l = 1,llm
     473
     474               pgcm = aps(l) + bps(l)*ps(i,j)  ! gcm pressure
     475               pgcm = log(pgcm)                ! for the vertical interpolation
     476               mmean(i,j,l) = 0.
     477
     478! 3.1 vertical interpolation
     479
     480               do n = 1,nspe
     481                  call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt)
     482                  vmrint(n) = vmr
     483                  iq = niq(n)
     484                  mmean(i,j,l) = mmean(i,j,l) + vmrint(n)*mmol(iq)
     485               end do
     486
     487! 3.2 attribute mixing ratio: - all layers or only thermosphere
     488!                             - with our without h2o
     489
     490               if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 1.e-3)) then
     491                  do n = 1,nspe
     492                     iq = niq(n)
     493                     if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
     494                        pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
     495                     end if
     496                  end do
     497               end if
     498
     499            end do
     500         end do
     501      end do
     502
     503      ! set surface values of chemistry tracers to zero
     504      if (flagthermo == 0) then
     505        ! NB: no problem for "surface water vapour" tracer which is always 0
     506        do n=1,nspe
     507          iq=niq(n)
     508          qsurf(1:ngridmx,iq)=0
    411509        enddo
    412510      endif
    413511
    414       ! build local chemistry tracer index array nqchem(:)
    415       ! as in the old days, water vapour is last and water ice,
    416       ! if present, is just before water vapour
    417       do iq=1,16 ! loop so as to avoid out-of-bounds on array
    418       if (iq==1) nqchem(i)=igcm_co2
    419       if (iq==2) nqchem(i)=igcm_co
    420       if (iq==3) nqchem(i)=igcm_o
    421       if (iq==4) nqchem(i)=igcm_o1d
    422       if (iq==5) nqchem(i)=igcm_o2
    423       if (iq==6) nqchem(i)=igcm_o3
    424       if (iq==7) nqchem(i)=igcm_h
    425       if (iq==8) nqchem(i)=igcm_h2
    426       if (iq==9) nqchem(i)=igcm_oh
    427       if (iq==10) nqchem(i)=igcm_ho2
    428       if (iq==11) nqchem(i)=igcm_h2o2
    429       if (iq==12) nqchem(i)=igcm_ch4
    430       if (iq==13) nqchem(i)=igcm_n2
    431       if (iq==14) nqchem(i)=igcm_ar
    432       if (iq==15) nqchem(i)=igcm_h2o_ice
    433       if (iq==16) nqchem(i)=igcm_h2o_vap
    434       enddo
    435 
    436 ! Load in chemistry data for initialization:
    437 c----------------------------------------------------------------------
    438 c Order of Major species in file:
    439 c
    440 c    1=CO2
    441 c    2=AR
    442 c    3=N2 
    443 c    4=O2 
    444 c    5=CO 
    445 c    6=O   
    446 c    7=H2
    447 c
    448 c Order of Minor species in files:
    449 c
    450 c    1=H 
    451 c    2=OH
    452 c    3=HO2
    453 c    4=H2O
    454 c    5=H2O2
    455 c    6=O1D
    456 c    7=O3
    457 c
    458 c----------------------------------------------------------------------
    459 c Major species:
    460         aa(igcm_co2) = 1
    461         aa(igcm_ar)  = 2
    462         aa(igcm_n2)  = 3
    463         aa(igcm_o2)  = 4
    464         aa(igcm_co)  = 5
    465         aa(igcm_o)   = 6
    466         aa(igcm_h2)  = 7
    467 c Minor species:
    468         aa(igcm_h)       = 8
    469         aa(igcm_oh)      = 9
    470         aa(igcm_ho2)     = 10
    471         aa(igcm_h2o_vap) = 11
    472         aa(igcm_h2o2)    = 12
    473         aa(igcm_o1d)     = 13
    474         aa(igcm_o3)      = 14
    475 c----------------------------------------------------------------------
    476       open(210, iostat=ierr,file=
    477      & trim(datafile)//'/atmosfera_LMD_may.dat')
    478       if (ierr.ne.0) then
    479         write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
    480         write(*,*)'(in aeronomars/inichim.F)'
    481         write(*,*)'It should be in :', trim(datafile),'/'
    482         write(*,*)'1) You can change this path in callphys.def with'
    483         write(*,*)'   datadir=/path/to/datafiles/'
    484         write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
    485         write(*,*)'   can be obtained online on:'
    486         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    487          STOP
    488       endif
    489       open(220, iostat=ierr,file=
    490      & trim(datafile)//'/atmosfera_LMD_min.dat')
    491       if (ierr.ne.0) then
    492         write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
    493         write(*,*)'(in aeronomars/inichim.F)'
    494         write(*,*)'It should be in :', trim(datafile),'/'
    495         write(*,*)'1) You can change this path in callphys.def with'
    496         write(*,*)'   datadir=/path/to/datafiles/'
    497         write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
    498         write(*,*)'   can be obtained online on:'
    499         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
    500          STOP
    501       endif
    502 
    503       read(210,*) tmp
    504       read(220,*) tmp
    505         do l=1,252
    506           read(210,112) tfile(l),zfile(l),nxf(l),(vmrf(l,n),n=1,7)
    507           zfile(l)=zfile(l)*100.              !pressure (Pa)
    508           sigfile(l)=zfile(l)/zfile(1)       
    509         enddo
    510         do l=1,252
    511           read(220,113)zfilemin(l),(vmrf(l,n),n=8,14)
    512           zfilemin(l)=zfilemin(l)*1000.       !height (m)
    513           do n=1,14
    514             nf(l,n)=vmrf(l,n)*nxf(l)
    515           enddo
    516         enddo
    517         close(210)
    518         close(220)
    519 
    520 c flag thermo set to 1 or 0 in newstart
    521 c inil=33 for initialisation of thermosphere only       
    522 c inil=1 for initialisation of all layers       
    523         if (thermo.eq.1) then
    524         inil=33
    525         else
    526         inil=1
    527         endif
    528 
    529 ! Initialization of tracers
    530 
    531       do i=1,iip1
    532        do j=1,jjp1
    533         do l=inil,llm
    534 
    535 c          zgcm=sig(l)
    536           zgcm=sig(l)*ps(i,j)
    537           densconc=0.
    538           nt=0.
    539 
    540           do n=1,14
    541             call intrplf(zgcm,vmrint,zfile,nf(1,n),252)
    542 c            call intrplf(zgcm,vmrint,sigfile,nf(1,n),252)
    543             ni(n)=vmrint
    544 
    545             densconc=densconc+ni(n)*mmolf(n)
    546             nt=nt+ni(n)
    547           enddo
    548 
    549           mmean=densconc/nt                     ! in amu
    550 
    551           if (nq.ne.nbqchem) then ! don't initialize water vapour
    552             do n=1,nq-1
    553             pq(i,j,l,nqchem(n))=
    554      &                 ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean
    555             if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n))
    556             enddo
    557             if (water .and. (nq.gt.nbqchem-2)) then
    558               pq(i,j,l,nqchem(nq)) = 0.
    559               if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq))
    560             else
    561               pq(i,j,l,nqchem(nq))=
    562      &              ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
    563               if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq))
    564             endif
    565           endif ! of if (nq.ne.nbqchem)
    566          
    567           if (nq.eq.nbqchem) then ! also initialize water vapour
    568             do n=1,nq-2
    569               pq(i,j,l,nqchem(n))=
    570      &                   ni(aa(nqchem(n)))/nt*mmol(nqchem(n))/mmean
    571               if(i.eq.iip1) pq(i,j,l,nqchem(n))=pq(1,j,l,nqchem(n))
    572             enddo
    573             if (water) then
    574               pq(i,j,l,nqchem(nq-1)) = 0.
    575               if(i.eq.iip1) then
    576                 pq(i,j,l,nqchem(nq-1))=pq(1,j,l,nqchem(nq-1))
    577               endif
    578               pq(i,j,l,nqchem(nq))=
    579      &                   ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
    580               if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq))
    581             else
    582               do n=nq-1,nq
    583                 pq(i,j,l,nqchem(nq))=
    584      &                   ni(aa(nqchem(nq)))/nt*mmol(nqchem(nq))/mmean
    585                 if(i.eq.iip1) pq(i,j,l,nqchem(nq))=pq(1,j,l,nqchem(nq))
    586               enddo
    587             endif ! of if (water)
    588           endif ! of if (nq.eq.nqmx)
    589         enddo     !nlayer
    590        enddo      !ngrid_j
    591       enddo       !ngrid_i
    592 
    593 cccccccccccccccccccccccccccccc     
    594 c  make sure that sum of q = 1     
    595 c dominent species is = 1 - sum(all other species)     
    596 cccccccccccccccccccccccccccccc     
    597 !      iqmax=nqchem_min
    598        iqmax=1
    599      
    600 !      if ((nqmx-nqchem_min).gt.10) then
    601       if (nbqchem.gt.10) then
    602        do l=1,llm
    603         do j=1,jjp1
    604           do i=1,iip1
    605 !           do iq=nqchem_min,nqmx
    606             do iiq=1,nbqchem
    607               iq=nqchem(iiq)
    608               if ( (pq(i,j,l,iq).gt.pq(i,j,l,iqmax)).and.
    609      .           (noms(iq).ne."ice") )  then
    610               iqmax=iq
    611               endif
    612             enddo
    613             pq(i,j,l,iqmax)=1.
    614             qtot(i,j,l)=0
    615 !           do iq=nqchem_min,nqmx
    616             do iiq=1,nbqchem
    617              iq=nqchem(iiq)
    618              if ( (iq.ne.iqmax).and.
    619      .           (noms(iq).ne."ice") ) then       
    620                pq(i,j,l,iqmax)=pq(i,j,l,iqmax)-pq(i,j,l,iq)       
    621              endif
    622             enddo !iq
    623 !            do iq=nqchem_min,nqmx
    624             do iiq=1,nbqchem
    625              iq=nqchem(iiq)
    626               if (noms(iq).ne."ice") then
    627                 qtot(i,j,l)=qtot(i,j,l)+pq(i,j,l,iq)
    628               endif
    629 c            if (i.eq.1.and.j.eq.1.and.l.eq.1) write(*,*) 'qtot(i,j,l)',
    630 c     $    qtot(i,j,l)
    631             enddo !iq
    632             if (i.eq.1.and.j.eq.1) write(*,*) 'inichim_newstart: ',
    633      $        'qtot(i,j,l)=',qtot(i,j,l)
    634           enddo !i   
    635          enddo !j   
    636        enddo !l 
    637       endif ! of if (nbqchem.gt.10)
    638 ccccccccccccccccccccccccccccccc
    639 
    640 66      format(i2,6(1x,e11.5))
    641 112     format(7x,f7.3,3x,e12.6,3x,e12.6,7(3x,e12.6))
    642 113     format(1x,f6.2,7(3x,e12.6))
     512
     513! 3.3 initialization of tracers not contained in the initialization files
     514
     515! methane : 10 ppbv
     516
     517      if (igcm_ch4 /= 0) then
     518         vmr = 10.e-9       
     519         do i = 1,iip1
     520            do j = 1,jjp1
     521               do l = 1,llm
     522                  pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l)
     523               end do
     524            end do
     525         end do
     526         ! set surface value to zero
     527         qsurf(1:ngridmx,igcm_ch4)=0
     528      end if
    643529
    644530      end

  • trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F

    r576 r618  
    6060      parameter (nesp = 16)      ! number of species in the chemistry
    6161c
    62       real ctimestep             ! chemistry timestep (s)
     62      real stimestep             ! standard timestep for the chemistry (s)
     63      real ctimestep             ! real timestep for the chemistry (s)
    6364      real rm(nlayermx,nesp)     ! species volume mixing ratio
    6465      real j(nlayermx,nd)        ! interpolated photolysis rates (s-1)
    … …  
    8586c
    8687cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    87 c     ctimestep : chemistry timestep (s)                         c
    88 c                 1/3 of physical timestep                       c
    89 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    90 c
    91       phychemrat = 3
     88c     stimestep  : standard timestep for the chemistry (s)       c
     89c     ctimestep  : real timestep for the chemistry (s)           c
     90c     phychemrat : ratio physical/chemical timestep              c
     91cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     92c
     93      stimestep = 600. ! standard value : 10 mn
     94c
     95      phychemrat = nint(ptimestep/stimestep)
    9296c
    9397      ctimestep = ptimestep/real(phychemrat)
    … …  
    11751179c
    11761180      integer l,iq
    1177       logical,save :: firstcall = .true.
    11781181     
    11791182c     tracer indexes in the chemistry:
    … …  
    11961199      integer,parameter :: i_ox   = 16
    11971200c     
    1198 c     first call initializations
    1199 c
    1200       if (firstcall) then
    1201 
    1202 c       identify the indexes of the tracers we need
    1203 
    1204         if (igcm_co2.eq.0) then
    1205           write(*,*) "concentrations: Error; no CO2 tracer !!!"
    1206           stop
    1207         endif
    1208         if (igcm_co.eq.0) then
    1209           write(*,*) "concentrations: Error; no CO tracer !!!"
    1210           stop
    1211         endif
    1212         if (igcm_o.eq.0) then
    1213           write(*,*) "concentrations: Error; no O tracer !!!"
    1214           stop
    1215         endif
    1216         if (igcm_o1d.eq.0) then
    1217           write(*,*) "concentrations: Error; no O1D tracer !!!"
    1218           stop
    1219         endif
    1220         if (igcm_o2.eq.0) then
    1221           write(*,*) "concentrations: Error; no O2 tracer !!!"
    1222           stop
    1223         endif
    1224         if (igcm_o3.eq.0) then
    1225           write(*,*) "concentrations: Error; no O3 tracer !!!"
    1226           stop
    1227         endif
    1228         if (igcm_h.eq.0) then
    1229           write(*,*) "concentrations: Error; no H tracer !!!"
    1230           stop
    1231         endif
    1232         if (igcm_h2.eq.0) then
    1233           write(*,*) "concentrations: Error; no H2 tracer !!!"
    1234           stop
    1235         endif
    1236         if (igcm_oh.eq.0) then
    1237           write(*,*) "concentrations: Error; no OH tracer !!!"
    1238           stop
    1239         endif
    1240         if (igcm_ho2.eq.0) then
    1241           write(*,*) "concentrations: Error; no HO2 tracer !!!"
    1242           stop
    1243         endif
    1244         if (igcm_h2o2.eq.0) then
    1245           write(*,*) "concentrations: Error; no H2O2 tracer !!!"
    1246           stop
    1247         endif
    1248         if (igcm_ch4.eq.0) then
    1249           write(*,*) "concentrations: Error; no CH4 tracer !!!"
    1250           stop
    1251         endif
    1252         if (igcm_n2.eq.0) then
    1253           write(*,*) "concentrations: Error; no N2 tracer !!!"
    1254           stop
    1255         endif
    1256         if (igcm_ar.eq.0) then
    1257           write(*,*) "concentrations: Error; no AR tracer !!!"
    1258           stop
    1259         endif
    1260         if (igcm_h2o_vap.eq.0) then
    1261           write(*,*) "concentrations: Error; no water vapor tracer !!!"
    1262           stop
    1263         endif
    1264         firstcall = .false.
    1265       end if ! of if (firstcall)
    1266 c
    12671201cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12681202c     initialise chemical species
    … …  
    12811215         rm(l,i_ho2)  = max(zycol(l, igcm_ho2),     1.e-30)
    12821216         rm(l,i_h2o2) = max(zycol(l, igcm_h2o2),    1.e-30)
    1283          rm(l,i_ch4)  = max(zycol(l, igcm_ch4),     1.e-30)
    12841217         rm(l,i_n2)   = max(zycol(l, igcm_n2),      1.e-30)
    12851218         rm(l,i_h2o)  = max(zycol(l, igcm_h2o_vap), 1.e-30)
    12861219      end do
     1220
     1221      if (igcm_ch4 .eq. 0) then
     1222         do l = 1,lswitch-1
     1223            rm(l,i_ch4) = 0.
     1224         end do
     1225      else
     1226         do l = 1,lswitch-1
     1227            rm(l,i_ch4) = max(zycol(l,igcm_ch4), 1.e-30)
     1228         end do
     1229      end if
    12871230c
    12881231cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    13701313         zycol(l, igcm_ho2)     = rm(l,i_ho2)
    13711314         zycol(l, igcm_h2o2)    = rm(l,i_h2o2)
    1372          zycol(l, igcm_ch4)     = rm(l,i_ch4)
    13731315         zycol(l, igcm_n2)      = rm(l,i_n2)
    13741316         zycol(l, igcm_h2o_vap) = rm(l,i_h2o)
    13751317      end do
     1318
     1319      if (igcm_ch4 .ne. 0) then
     1320         do l = 1,lswitch-1
     1321            zycol(l,igcm_ch4) = rm(l,i_ch4)
     1322         end do
     1323      end if
    13761324c
    13771325      return
    … …  
    20001948      return
    20011949      end
    2002 
Note: See TracChangeset for help on using the changeset viewer.