Changeset 1125


Ignore:
Timestamp:
Dec 11, 2013, 5:43:59 PM (11 years ago)
Author:
flefevre
Message:

calchim.F90 et photochemistry.F:

Suppression des references a nlayermx et nqmx.
=> Passage de nlayer et nq en arguments de subroutines.

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

Legend:

Unmodified
Added
Removed

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

    r1047 r1125  
    1313                            igcm_noplus, igcm_n2plus, igcm_hplus,         &
    1414                            igcm_hco2plus, igcm_elec, mmol
     15
    1516      use conc_mod, only: mmean ! mean molecular mass of the atmosphere
    1617
    … …  
    3334!    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
    3435!    update 16/03/2012 optimization (Franck Lefevre)
     36!    update 11/12/2013 optimization (Franck Lefevre)
    3537!
    3638!   Arguments:
    … …  
    7476!     input:
    7577
    76       integer,intent(in) :: ngrid ! number of atmospheric columns
    77       integer,intent(in) :: nlayer ! number of atmospheric layers
    78       integer,intent(in) :: nq ! number of tracers
     78      integer,intent(in) :: ngrid    ! number of atmospheric columns
     79      integer,intent(in) :: nlayer   ! number of atmospheric layers
     80      integer,intent(in) :: nq       ! number of tracers
    7981      real :: ptimestep
    8082      real :: pplay(ngrid,nlayer)    ! pressure at the middle of the layers
    … …  
    8890      real :: pv(ngrid,nlayer)       ! v component of the wind (m.s-1)
    8991      real :: pdv(ngrid,nlayer)      ! v component tendency
    90       real :: dist_sol                   ! distance of the sun (AU)
    91       real :: mu0(ngrid)               ! cos of solar zenith angle (=1 when sun at zenith)
    92       real :: pq(ngrid,nlayer,nq)  ! tracers mass mixing ratio
    93       real :: pdq(ngrid,nlayer,nq) ! previous tendencies
    94       real :: zday                       ! date (time since Ls=0, in martian days)
    95       real :: tauref(ngrid)            ! optical depth at 7 hPa
    96       real :: co2ice(ngrid)            ! co2 ice surface layer (kg.m-2)
     92      real :: dist_sol               ! distance of the sun (AU)
     93      real :: mu0(ngrid)             ! cos of solar zenith angle (=1 when sun at zenith)
     94      real :: pq(ngrid,nlayer,nq)    ! tracers mass mixing ratio
     95      real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
     96      real :: zday                   ! date (time since Ls=0, in martian days)
     97      real :: tauref(ngrid)          ! optical depth at 7 hPa
     98      real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
    9799      real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
    98100      real :: surfice(ngrid,nlayer)  !  ice surface area (m2/m3)
    … …  
    100102!     output:
    101103
    102       real :: dqchim(ngrid,nlayer,nq) ! tendencies on pq due to chemistry
     104      real :: dqchim(ngrid,nlayer,nq)   ! tendencies on pq due to chemistry
    103105      real :: dqschim(ngrid,nq)         ! tendencies on qsurf
    104       real :: dqcloud(ngrid,nlayer,nq)! tendencies on pq due to condensation
     106      real :: dqcloud(ngrid,nlayer,nq)  ! tendencies on pq due to condensation
    105107      real :: dqscloud(ngrid,nq)        ! tendencies on qsurf
    106108
    … …  
    109111      integer,save :: nbq                   ! number of tracers used in the chemistry
    110112      integer,allocatable,save :: niq(:)    ! array storing the indexes of the tracers
    111       integer :: iloc(1)            ! index of major species
     113      integer :: iloc(1)                    ! index of major species
    112114      integer :: ig,l,i,iq,iqmax
    113115      integer :: foundswitch, lswitch
    … …  
    147149
    148150      real    :: latvl1, lonvl1
    149       real    :: zq(ngrid,nlayer,nq) ! pq+pdq*ptimestep before chemistry
    150                                            ! new mole fraction after
     151      real    :: zq(ngrid,nlayer,nq)   ! pq+pdq*ptimestep before chemistry
     152                                       ! new mole fraction after
    151153      real    :: zt(ngrid,nlayer)      ! temperature
    152154      real    :: zu(ngrid,nlayer)      ! u component of the wind
    153155      real    :: zv(ngrid,nlayer)      ! v component of the wind
    154       real    :: taucol                    ! optical depth at 7 hPa
     156      real    :: taucol                ! optical depth at 7 hPa
    155157
    156158      logical,save :: firstcall = .true.
    157       logical,save :: depos = .false.      ! switch for dry deposition
     159      logical,save :: depos = .false.  ! switch for dry deposition
    158160
    159161!     for each column of atmosphere:
    160162
    161       real :: zpress(nlayer)       !  Pressure (mbar)
    162       real :: zdens(nlayer)        !  Density  (cm-3)
    163       real :: ztemp(nlayer)        !  Temperature (K)
    164       real :: zlocal(nlayer)       !  Altitude (km)
    165       real :: zycol(nlayer,nq)   ! Composition (mole fractions)
    166       real :: szacol                 ! Solar zenith angle
    167       real :: surfice1d(nlayer)    !  Ice surface area (cm2/cm3)
    168       real :: surfdust1d(nlayer)   !  Dust surface area (cm2/cm3)
    169       real :: jo3(nlayer)          !  Photodissociation rate O3->O1D (s-1)
     163      real :: zpress(nlayer)       ! Pressure (mbar)
     164      real :: zdens(nlayer)        ! Density  (cm-3)
     165      real :: ztemp(nlayer)        ! Temperature (K)
     166      real :: zlocal(nlayer)       ! Altitude (km)
     167      real :: zycol(nlayer,nq)     ! Composition (mole fractions)
     168      real :: szacol               ! Solar zenith angle
     169      real :: surfice1d(nlayer)    ! Ice surface area (cm2/cm3)
     170      real :: surfdust1d(nlayer)   ! Dust surface area (cm2/cm3)
     171      real :: jo3(nlayer)          ! Photodissociation rate O3->O1D (s-1)
    170172
    171173!     for output:
    172174
    173       logical :: output                 ! to issue calls to writediagfi and stats
     175      logical :: output             ! to issue calls to writediagfi and stats
    174176      parameter (output = .true.)
    175177      real :: jo3_3d(ngrid,nlayer)  ! Photodissociation rate O3->O1D (s-1)
    … …  
    645647
    646648         if (photochem) then
    647             call photochemistry(lswitch,zycol,szacol,ptimestep,    &
     649            call photochemistry(nlayer,nq,                         &
     650                                lswitch,zycol,szacol,ptimestep,    &
    648651                                zpress,ztemp,zdens,dist_sol,       &
    649652                                surfdust1d,surfice1d,jo3,taucol)

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

    r1036 r1125  
    66c     ------
    77c
    8 c     Version: 17/03/2011
     8c     Version: 11/12/2013
    99c
    1010c*****************************************************************
    11 c
    12       subroutine photochemistry(lswitch, zycol, sza, ptimestep,
     11
     12      subroutine photochemistry(nlayer, nq,
     13     $                  lswitch, zycol, sza, ptimestep,
    1314     $                  press,temp, dens, dist_sol, surfdust1d,
    1415     $                  surfice1d, jo3, tau)
    15 c
    16       use tracer_mod, only: nqmx
     16
    1717      implicit none
    18 c
    19 #include "dimensions.h"
    20 #include "dimphys.h"
     18
     19!#include "dimensions.h"
     20!#include "dimphys.h"
    2121#include "chimiedata.h"
    2222#include "callkeys.h"
    2323!#include "tracer.h"
    24 c
     24
     25cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     26c     inputs:
     27cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     28
     29      integer, intent(in) :: nlayer ! number of atmospheric layers
     30      integer, intent(in) :: nq     ! number of tracers
     31     
     32      real :: sza                   ! solar zenith angle (deg)
     33      real :: ptimestep             ! physics timestep (s)
     34      real :: press(nlayer)         ! pressure (hpa)
     35      real :: temp(nlayer)          ! temperature (k)
     36      real :: dens(nlayer)          ! density (cm-3)
     37      real :: dist_sol              ! sun distance (au)
     38      real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
     39      real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
     40      real :: tau                   ! optical depth at 7 hpa
     41
    2542cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    2643c     input/output:
    2744cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    28 c     
    29       real zycol(nlayermx,nqmx)  ! chemical species volume mixing ratio
    30 c
    31 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    32 c     inputs:
    33 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    34 c     
    35       real sza                   ! solar zenith angle (deg)
    36       real ptimestep             ! physics timestep (s)
    37       real press(nlayermx)       ! pressure (hpa)
    38       real temp(nlayermx)        ! temperature (k)
    39       real dens(nlayermx)        ! density (cm-3)
    40       real dist_sol              ! sun distance (au)
    41       real surfdust1d(nlayermx)  ! dust surface area (cm2/cm3)
    42       real surfice1d(nlayermx)   ! ice surface area (cm2/cm3)
    43       real tau                   ! optical depth at 7 hpa
    44 c
     45     
     46      real :: zycol(nlayer,nq)      ! chemical species volume mixing ratio
     47
    4548cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    4649c     output:
    4750cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    48 c     
    49       real jo3(nlayermx)  ! photodissociation rate o3 -> o1d
    50 c
     51     
     52      real :: jo3(nlayer)  ! photodissociation rate o3 -> o1d
     53
    5154cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    5255c     local:
    5356cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    54 c
    55       integer phychemrat         ! ratio physics timestep/chemistry timestep
    56       integer istep
    57       integer i_co2,i_o3,j_o3_o1d,lev
    58       integer nesp
    59       integer lswitch
    60 c
    61       parameter (nesp = 16)      ! number of species in the chemistry
    62 c
    63       real stimestep             ! standard timestep for the chemistry (s)
    64       real ctimestep             ! real timestep for the chemistry (s)
    65       real rm(nlayermx,nesp)     ! species volume mixing ratio
    66       real j(nlayermx,nd)        ! interpolated photolysis rates (s-1)
    67       real rmco2(nlayermx)       ! co2 mixing ratio
    68       real rmo3(nlayermx)        ! ozone mixing ratio
    69 c
     57
     58      integer, parameter :: nesp = 16 ! number of species in the chemistry
     59
     60      integer :: phychemrat       ! ratio physics timestep/chemistry timestep
     61      integer :: istep
     62      integer :: i_co2,i_o3,j_o3_o1d,lev
     63      integer :: lswitch
     64
     65      real :: stimestep           ! standard timestep for the chemistry (s)
     66      real :: ctimestep           ! real timestep for the chemistry (s)
     67      real :: rm(nlayer,nesp)     ! species volume mixing ratio
     68      real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
     69      real :: rmco2(nlayer)       ! co2 mixing ratio
     70      real :: rmo3(nlayer)        ! ozone mixing ratio
     71
    7072c     reaction rates
    71 c
    72       real a001(nlayermx), a002(nlayermx), a003(nlayermx)
    73       real b001(nlayermx), b002(nlayermx), b003(nlayermx),
    74      $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
    75      $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
    76       real c001(nlayermx), c002(nlayermx), c003(nlayermx),
    77      $     c004(nlayermx), c005(nlayermx), c006(nlayermx),
    78      $     c007(nlayermx), c008(nlayermx), c009(nlayermx),
    79      $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
    80      $     c013(nlayermx), c014(nlayermx), c015(nlayermx),
    81      $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
    82       real d001(nlayermx), d002(nlayermx), d003(nlayermx)
    83       real e001(nlayermx), e002(nlayermx), e003(nlayermx)
    84       real h001(nlayermx), h002(nlayermx), h003(nlayermx),
    85      $     h004(nlayermx), h005(nlayermx)
    86       real t001(nlayermx), t002(nlayermx), t003(nlayermx)
    87 c
     73
     74      real :: a001(nlayer), a002(nlayer), a003(nlayer)
     75      real :: b001(nlayer), b002(nlayer), b003(nlayer),
     76     $        b004(nlayer), b005(nlayer), b006(nlayer),
     77     $        b007(nlayer), b008(nlayer), b009(nlayer)
     78      real :: c001(nlayer), c002(nlayer), c003(nlayer),
     79     $        c004(nlayer), c005(nlayer), c006(nlayer),
     80     $        c007(nlayer), c008(nlayer), c009(nlayer),
     81     $        c010(nlayer), c011(nlayer), c012(nlayer),
     82     $        c013(nlayer), c014(nlayer), c015(nlayer),
     83     $        c016(nlayer), c017(nlayer), c018(nlayer)
     84      real :: d001(nlayer), d002(nlayer), d003(nlayer)
     85      real :: e001(nlayer), e002(nlayer), e003(nlayer)
     86      real :: h001(nlayer), h002(nlayer), h003(nlayer),
     87     $        h004(nlayer), h005(nlayer)
     88      real :: t001(nlayer), t002(nlayer), t003(nlayer)
     89
    8890cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    8991c     stimestep  : standard timestep for the chemistry (s)       c
    … …  
    9193c     phychemrat : ratio physical/chemical timestep              c
    9294cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    93 c
     95
    9496      stimestep = 600. ! standard value : 10 mn
    95 c
     97
    9698      phychemrat = nint(ptimestep/stimestep)
    97 c
     99
    98100      ctimestep = ptimestep/real(phychemrat)
    99 c
     101
    100102cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    101103c     initialisation of chemical species and families            c
    102104cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    103 c
    104       call gcmtochim(zycol, lswitch, nesp, rm)
    105 c
     105
     106      call gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
     107
    106108cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    107109c     compute reaction rates                                     c
    108110cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    109 c
    110       call chemrates(lswitch, dens, press, temp,
     111
     112      call chemrates(nlayer,
     113     $               lswitch, dens, press, temp,
    111114     $               surfdust1d, surfice1d,
    112115     $               a001, a002, a003,
    … …  
    120123     $               h001, h002, h003, h004, h005,
    121124     $               t001, t002, t003, tau)
    122 c
     125
    123126cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    124127c     interpolation of photolysis rates in the lookup table      c
    125128cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    126 c
     129
    127130      i_co2 = 1
    128131      i_o3  = 6
    129 c
     132
    130133      do lev = 1,lswitch-1
    131134         rmco2(lev) = rm(lev,i_co2)
    132135         rmo3(lev)  = rm(lev, i_o3)
    133136      end do
    134 c
    135       call phot(lswitch, press, temp, sza, tau, dist_sol,
     137
     138      call phot(nlayer, lswitch, press, temp, sza, tau, dist_sol,
    136139     $          rmco2, rmo3, j)
    137 c
     140
    138141      j_o3_o1d = 5
    139 c
     142
    140143      do lev = 1,lswitch-1
    141144         jo3(lev) = j(lev,j_o3_o1d)
    142145      end do
    143 c
     146
    144147cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    145148c     loop over time within the physical timestep                c
    146149cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    147 c
     150
    148151      do istep = 1,phychemrat
    149          call chimie(lswitch,nesp, rm, j, dens, ctimestep,
     152         call chimie(nlayer, lswitch, nesp, rm, j, dens, ctimestep,
    150153     $               press, temp, sza, dist_sol,
    151154     $               a001, a002, a003,
    … …  
    160163     $               t001, t002, t003)
    161164      end do
    162 c
     165
    163166cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    164167c     save chemical species for the gcm                          c
    165168cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    166 c
    167       call chimtogcm(zycol, lswitch, nesp, rm)
    168 c
     169
     170      call chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
     171
    169172      return
    170173      end
    171 c
     174
    172175c*****************************************************************
    173 c
    174       subroutine chimie(lswitch, nesp, rm, j, dens, dt,
     176
     177      subroutine chimie(nlayer,
     178     $                  lswitch, nesp, rm, j, dens, dt,
    175179     $                  press, t, sza, dist_sol,
    176180     $                  a001, a002, a003,
    … …  
    184188     $                  h001, h002, h003, h004, h005,
    185189     $                  t001, t002, t003)
    186 c
     190
    187191c*****************************************************************
    188 c
     192
    189193      implicit none
    190 c
    191 #include "dimensions.h"
    192 #include "dimphys.h"
     194
     195!#include "dimensions.h"
     196!#include "dimphys.h"
    193197#include "chimiedata.h"
    194198#include "callkeys.h"
    195 c
     199
     200cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     201c     inputs:
     202cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     203
     204      integer, intent(in) :: nlayer ! number of atmospheric layers
     205      integer :: lswitch            ! interface level between chemistries
     206      integer :: nesp               ! number of species
     207
     208      real :: dens(nlayer)          ! density (cm-3)
     209      real :: dt                    ! chemistry timestep (s)
     210      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1)
     211      real :: press(nlayer)         ! pressure (hpa)
     212      real :: t(nlayer)             ! temperature (k)
     213      real :: sza                   ! solar zenith angle (deg)
     214      real :: dist_sol              ! sun distance (au)
     215
     216c     reaction rates
     217
     218      real :: a001(nlayer), a002(nlayer), a003(nlayer)
     219      real :: b001(nlayer), b002(nlayer), b003(nlayer),
     220     $        b004(nlayer), b005(nlayer), b006(nlayer),
     221     $        b007(nlayer), b008(nlayer), b009(nlayer)
     222      real :: c001(nlayer), c002(nlayer), c003(nlayer),
     223     $        c004(nlayer), c005(nlayer), c006(nlayer),
     224     $        c007(nlayer), c008(nlayer), c009(nlayer),
     225     $        c010(nlayer), c011(nlayer), c012(nlayer),
     226     $        c013(nlayer), c014(nlayer), c015(nlayer),
     227     $        c016(nlayer), c017(nlayer), c018(nlayer)
     228      real :: d001(nlayer), d002(nlayer), d003(nlayer)
     229      real :: e001(nlayer), e002(nlayer), e003(nlayer)
     230      real :: h001(nlayer), h002(nlayer), h003(nlayer),
     231     $        h004(nlayer), h005(nlayer)
     232      real :: t001(nlayer), t002(nlayer), t003(nlayer)
     233
    196234cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    197235c     input/output:
    198236cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    199 c     
    200       integer lswitch         ! interface level between chemistries
    201       integer nesp            ! number of species
    202 c
    203       real rm(nlayermx,nesp)  ! volume mixing ratios
    204 c
    205 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    206 c     inputs:
    207 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    208 c
    209       real dens(nlayermx)     ! density (cm-3)
    210       real dt                 ! chemistry timestep (s)
    211       real j(nlayermx,nd)     ! interpolated photolysis rates (s-1)
    212       real press(nlayermx)    ! pressure (hpa)
    213       real t(nlayermx)        ! temperature (k)
    214       real sza                ! solar zenith angle (deg)
    215       real dist_sol           ! sun distance (au)
    216 c
    217 c     reaction rates
    218 c
    219       real a001(nlayermx), a002(nlayermx), a003(nlayermx)
    220       real b001(nlayermx), b002(nlayermx), b003(nlayermx),
    221      $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
    222      $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
    223       real c001(nlayermx), c002(nlayermx), c003(nlayermx),
    224      $     c004(nlayermx), c005(nlayermx), c006(nlayermx),
    225      $     c007(nlayermx), c008(nlayermx), c009(nlayermx),
    226      $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
    227      $     c013(nlayermx), c014(nlayermx), c015(nlayermx),
    228      $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
    229       real d001(nlayermx), d002(nlayermx), d003(nlayermx)
    230       real e001(nlayermx), e002(nlayermx), e003(nlayermx)
    231       real h001(nlayermx), h002(nlayermx), h003(nlayermx),
    232      $     h004(nlayermx), h005(nlayermx)
    233       real t001(nlayermx), t002(nlayermx), t003(nlayermx)
    234 c
     237
     238      real :: rm(nlayer,nesp)  ! volume mixing ratios
     239
    235240cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    236241c     local:
    237242cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    238 c
    239       real hetero_ice, hetero_dust
    240 c
    241       integer iesp, iter, l, niter
    242       integer i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
    243      $        i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
    244       integer j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
    245      $        j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
    246      $        j_ho2, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
    247      $        j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
    248      $        j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
    249      $        j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
    250      $        j_hno3, j_hno4
    251 c
    252       parameter (hetero_ice  = 1.)   ! switch for het. chem. on ice clouds
    253       parameter (hetero_dust = 0.)   ! switch for het. chem. on dust
    254                                      ! hetero_dust = 0. advised for the moment
    255 c
    256       parameter (niter = 5)          ! iterations in the chemical scheme
    257 c
    258       real cc0(nlayermx,nesp)        ! initial number density (cm-3)
    259       real cc(nlayermx,nesp)         ! final number density (cm-3)
    260       real nox(nlayermx)             ! nox number density (cm-3)
    261       real no(nlayermx)              ! no  number density (cm-3)
    262       real no2(nlayermx)             ! no2 number density (cm-3)
    263       real production(nlayermx,nesp) ! production rate
    264       real loss(nlayermx,nesp)       ! loss rate
    265 c
    266       real ro_o3, rh_ho2, roh_ho2, rno2_no
    267 c
     243
     244      real, parameter :: hetero_ice  = 1. ! switch for het. chem. on ice clouds
     245      real, parameter :: hetero_dust = 0. ! switch for het. chem. on dust
     246
     247      integer :: iesp, iter, l
     248      integer :: i_co2, i_co, i_o2, i_h2, i_h2o, i_h2o2, i_hox, i_ox,
     249     $           i_o1d, i_o, i_o3, i_h, i_oh, i_ho2, i_ch4, i_n2
     250      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d,
     251     $           j_o3_o1d, j_o3_o, j_h2o, j_hdo, j_h2o2,
     252     $           j_ho2, j_no2, j_ch4_ch3_h, j_ch4_1ch2_h2,
     253     $           j_ch4_3ch2_h_h, j_ch4_ch_h2_h, j_ch3o2h,
     254     $           j_ch2o_co, j_ch2o_hco, j_ch3oh, j_c2h6, j_hcl,
     255     $           j_hocl, j_clo, j_so2, j_so, j_h2s, j_so3,
     256     $           j_hno3, j_hno4
     257
     258      integer, parameter :: niter = 5        ! iterations in the chemical scheme
     259
     260      real :: cc0(nlayer,nesp)        ! initial number density (cm-3)
     261      real :: cc(nlayer,nesp)         ! final number density (cm-3)
     262      real :: nox(nlayer)             ! nox number density (cm-3)
     263      real :: no(nlayer)              ! no  number density (cm-3)
     264      real :: no2(nlayer)             ! no2 number density (cm-3)
     265      real :: production(nlayer,nesp) ! production rate
     266      real :: loss(nlayer,nesp)       ! loss rate
     267
     268      real :: ro_o3, rh_ho2, roh_ho2, rno2_no
     269
    268270cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    269271c     tracer numbering in the chemistry
    270272cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    271 c
     273
    272274      i_co2    =  1
    273275      i_co     =  2
    … …  
    286288      i_hox    = 15
    287289      i_ox     = 16
    288 c
     290
    289291cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    290292c     numbering of photolysis rates
    291293cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    292 c
     294
    293295      j_o2_o         =  1      ! o2 + hv     -> o + o
    294296      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
    … …  
    320322      j_hno3         =  28     ! hno3 + hv   -> oh + no2
    321323      j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
    322 c
     324
    323325cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    324326c     volume mixing ratio -> density conversion
    325327cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    326 c
     328
    327329      do iesp = 1,nesp
    328330         do l = 1,lswitch-1
    … …  
    331333         end do
    332334      end do
    333 c
     335
    334336cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    335337c     co2 and nox number densities (cm-3)   
    … …  
    337339c     nox mixing ratio: 6.e-10 (yung and demore, 1999)
    338340cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    339 c
     341
    340342      do l = 1,lswitch-1
    341343         nox(l) = 6.e-10*dens(l)
    342344      end do
    343 c
     345
    344346cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    345347c     loop over iterations in the chemical scheme
    346348cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    347 c
     349
    348350      do iter = 1,niter
    349 c
     351
    350352cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    351353c        nox species
    … …  
    353355c        no2/no
    354356cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    355 c
     357
    356358      do l = 1,lswitch-1
    357 c
     359
    358360         rno2_no = (d002(l)*cc(l,i_o3) + d003(l)*cc(l,i_ho2))
    359361     $            /(j(l,j_no2) +
    360362     $              d001(l)*max(cc(l,i_o),1.e-30*dens(l)))
    361 c
     363
    362364cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    363365c        no
    364366cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    365 c
     367
    366368         no(l) = nox(l)/(1. + rno2_no)
    367 c
     369
    368370cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    369371c        no2
    370372cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    371 c
     373
    372374         no2(l) = no(l)*rno2_no
    373 c
     375
    374376      end do
    375 c
     377
    376378cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    377379c        hox species
    … …  
    381383c        h/ho2
    382384cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    383 c
     385
    384386      do l = 1,lswitch-1
    385 c
     387
    386388         rh_ho2 = (c001(l)*cc(l,i_o)
    387389     $           + c004(l)*cc(l,i_h)
    … …  
    400402     $            /max(cc(l,i_h),dens(l)*1.e-30)   ! ajout 20090401
    401403     $           )
    402 c
     404
    403405cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    404406c        oh/ho2
    405407cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    406 c
     408
    407409         roh_ho2 = (c001(l)*cc(l,i_o)
    408410     $            + c003(l)*cc(l,i_o3)*rh_ho2
    … …  
    428430     $            + e001(l)*cc(l,i_co)
    429431     $            + h002(l)*hetero_ice)
    430 c
     432
    431433cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    432434c        h
    433435cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    434 c
     436
    435437         cc(l,i_h) = cc(l,i_hox)
    436438     $                 /(1. + (1. + roh_ho2)/rh_ho2)
    437 c
     439
    438440cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    439441c        ho2
    440442cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    441 c
     443
    442444         cc(l,i_ho2) = cc(l,i_h)/rh_ho2
    443 c
     445
    444446cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    445447c        oh
    446448cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    447 c
     449
    448450         cc(l,i_oh) = cc(l,i_ho2)*roh_ho2
    449 c
     451
    450452      end do
    451 c
     453
    452454cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    453455c        ox species
    … …  
    462464c        - continuity equation for o
    463465cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    464 c       
     466       
    465467      if (sza .le. 95.) then
    466 c
     468
    467469         do l = 1,lswitch-1
    468 c
     470
    469471cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    470472c        o(1d)
    471473cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    472 c
     474
    473475            cc(l,i_o1d) = (j(l,j_co2_o1d)*cc(l,i_co2)
    474476     $                   + j(l,j_o2_o1d)*cc(l,i_o2)
    … …  
    480482     $                   + b005(l)*cc(l,i_o3)
    481483     $                   + b006(l)*cc(l,i_o3))
    482 c
     484
    483485cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    484486c        o/o3
    485487cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    486 c
     488
    487489            ro_o3 = (j(l,j_o3_o1d) + j(l,j_o3_o)
    488490     $              + a003(l)*cc(l,i_o)
    … …  
    492494     $              )
    493495     $              /(a001(l)*cc(l,i_o2))
    494 c
     496
    495497cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    496498c        o3
    497499cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    498 c
     500
    499501            cc(l,i_o3) = cc(l,i_ox)/(1. + ro_o3)
    500 c
     502
    501503cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    502504c        o
    503505cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    504 c
     506
    505507            cc(l,i_o) = cc(l,i_o3)*ro_o3
    506 c
     508
    507509cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    508510c        ox = o + o3
    509511cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    510 c
     512
    511513            production(l,i_ox) =
    512514     $                   + j(l,j_co2_o)*cc(l,i_co2)
    … …  
    518520     $                   + c013(l)*cc(l,i_oh)*cc(l,i_oh)
    519521     $                   + d003(l)*cc(l,i_ho2)*no(l)
    520 c
     522
    521523            loss(l,i_ox) = 2.*a002(l)*cc(l,i_o)*cc(l,i_o)
    522524     $                   + 2.*a003(l)*cc(l,i_o)*cc(l,i_o3)
    … …  
    529531     $                   + d001(l)*cc(l,i_o)*no2(l)
    530532     $                   + e002(l)*cc(l,i_o)*cc(l,i_co)
    531 c
     533
    532534            loss(l,i_ox) = loss(l,i_ox)/cc(l,i_ox)
    533 c
     535
    534536         end do
    535 c
     537
    536538      else
    537 c
     539
    538540         do l = 1,lswitch-1
    539 c
     541
    540542cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    541543c        o(1d)
    542544cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    543 c
     545
    544546            cc(l,i_o1d) = 0.
    545 c
     547
    546548cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    547549c        o3
    548550cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    549 c
     551
    550552            production(l,i_o3) = a001(l)*cc(l,i_o2)*cc(l,i_o)
    551 c
     553
    552554            loss(l,i_o3) = a003(l)*cc(l,i_o)
    553555     $                   + c003(l)*cc(l,i_h)
    554556     $                   + c014(l)*cc(l,i_oh)
    555557     $                   + c015(l)*cc(l,i_ho2)
    556 c
    557558
    558559cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    559560c        o
    560561cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    561 c
     562
    562563            production(l,i_o) = c006(l)*cc(l,i_h)*cc(l,i_ho2)
    563564     $                        + c013(l)*cc(l,i_oh)*cc(l,i_oh)
    564 c
     565
    565566            loss(l,i_o)  =  a001(l)*cc(l,i_o2)
    566567     $                   + 2.*a002(l)*cc(l,i_o)
    … …  
    570571     $                   + c012(l)*cc(l,i_h2o2)
    571572     $                   + e002(l)*cc(l,i_co)
    572 c
     573
    573574         end do
    574575      end if
    575 c
     576
    576577cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    577578c     other species
    578579cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    579 c
     580
    580581      do l = 1,lswitch-1
    581 c
     582
    582583cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    583584c        co2
    584585cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    585 c
     586
    586587         production(l,i_co2) = e001(l)*cc(l,i_oh)*cc(l,i_co)
    587588     $                       + e002(l)*cc(l,i_o)*cc(l,i_co)
    588589     $                       + t002(l)*cc(l,i_ch4)*16./44. ! ajout 20090401
    589590     $                       + t003(l)*cc(l,i_co2)*8./44.  ! ajout 20090401
    590 c
     591
    591592         loss(l,i_co2) = j(l,j_co2_o)
    592593     $                 + j(l,j_co2_o1d)
    593594     $                 + t003(l)
    594 c
     595
    595596cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    596597c        co
    597598cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    598 c
     599
    599600         production(l,i_co) = j(l,j_co2_o)*cc(l,i_co2)
    600601     $                      + j(l,j_co2_o1d)*cc(l,i_co2)
    601602     $                      + t003(l)*cc(l,i_co2)
    602 c
     603
    603604         loss(l,i_co) = e001(l)*cc(l,i_oh)
    604605     $                + e002(l)*cc(l,i_o)
    605 c
     606
    606607cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    607608c        o2
    608609cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    609 c
     610
    610611         production(l,i_o2) =
    611612     $                  j(l,j_o3_o)*cc(l,i_o3)
    … …  
    625626     $                + c016(l)*cc(l,i_ho2)*cc(l,i_ho2)
    626627     $                + d001(l)*cc(l,i_o)*no2(l)
    627 c
     628
    628629         loss(l,i_o2) = j(l,j_o2_o)
    629630     $                + j(l,j_o2_o1d)
    630631     $                + a001(l)*cc(l,i_o)
    631632     $                + c011(l)*cc(l,i_h)
    632 c
     633
    633634cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    634635c        h2
    635636cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    636 c
     637
    637638         production(l,i_h2) = c005(l)*cc(l,i_h)*cc(l,i_ho2)
    638639     $                      + c018(l)*cc(l,i_h)*cc(l,i_h)
    639 c
     640
    640641         loss(l,i_h2) = b003(l)*cc(l,i_o1d)
    641642     $                + c010(l)*cc(l,i_oh)
    642 c
     643
    643644cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    644645c        h2o
    645646cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    646 c
     647
    647648         production(l,i_h2o) =
    648649     $                   c006(l)*cc(l,i_h)*cc(l,i_ho2)
    … …  
    652653     $                 + c013(l)*cc(l,i_oh)*cc(l,i_oh)
    653654     $                 + h004(l)*cc(l,i_h2o2)*hetero_ice
    654 c
     655
    655656         loss(l,i_h2o) = j(l,j_h2o)
    656657     $                 + b002(l)*cc(l,i_o1d)
    657658     $                 + t001(l)
    658 c
     659
    659660cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    660661c        h2o2
    661662cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    662 c
     663
    663664         production(l,i_h2o2) =
    664665     $                    c008(l)*cc(l,i_ho2)*cc(l,i_ho2)
    … …  
    667668c    $                  + 0.5*h001(l)*cc(l,i_ho2)*hetero_ice
    668669c    $                  + 0.5*h002(l)*cc(l,i_oh)*hetero_ice
    669 c
     670
    670671         loss(l,i_h2o2) = j(l,j_h2o2)
    671672     $                  + c009(l)*cc(l,i_oh)
    … …  
    673674     $                  + h004(l)*hetero_ice
    674675     $                  + h005(l)*hetero_dust
    675 c
     676
    676677cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    677678c        hox = h + oh + ho2
    678679cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    679 c
     680
    680681         production(l,i_hox) =
    681682     $                   2.*j(l,j_h2o)*cc(l,i_h2o)
    … …  
    685686     $                 + 2.*c012(l)*cc(l,i_o)*cc(l,i_h2o2)
    686687     $                 + 2.*t001(l)*cc(l,i_h2o)
    687 c
     688
    688689         loss(l,i_hox) = 2.*c005(l)*cc(l,i_h)*cc(l,i_ho2)
    689690     $                 + 2.*c006(l)*cc(l,i_h)*cc(l,i_ho2)
    … …  
    697698     $                 + h002(l)*cc(l,i_oh)*hetero_ice
    698699     $                 + h003(l)*cc(l,i_ho2)*hetero_dust
    699 c
     700
    700701         loss(l,i_hox) = loss(l,i_hox)/cc(l,i_hox)
    701 c
     702
    702703cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    703704c        ch4
    704705cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    705 c
     706
    706707         production(l,i_ch4) = 0.
    707 c
     708
    708709         loss(l,i_ch4) = j(l,j_ch4_ch3_h)
    709710     $                 + j(l,j_ch4_1ch2_h2)
    … …  
    714715     $                 + b009(l)*cc(l,i_o1d)
    715716     $                 + e003(l)*cc(l,i_oh)
    716 c
     717
    717718      end do
    718 c
     719
    719720cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    720721c     update number densities
    721722cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    722 c
     723
    723724c     long-lived species
    724 c
     725
    725726      do l = 1,lswitch-1
    726727         cc(l,i_co2) = (cc0(l,i_co2) + production(l,i_co2)*dt)
    … …  
    741742     $              /(1. + loss(l,i_ch4)*dt)
    742743      end do
    743 c
     744
    744745c     ox species
    745 c
     746
    746747      if (sza .le. 95.) then
    747748         do l = 1,lswitch-1
    … …  
    758759         end do
    759760      end if
    760 c
     761
    761762cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    762763c     end of loop over iterations
    763764cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    764 c
     765
    765766      end do
    766 c
     767
    767768cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    768769c     density -> volume mixing ratio conversion
    769770cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    770 c
     771
    771772      do iesp = 1,nesp
    772773         do l = 1,lswitch-1
    … …  
    774775         end do
    775776      end do
    776 c
     777
    777778      return
    778779      end
    779 c
     780
    780781c*****************************************************************
    781 c
    782       subroutine phot(lswitch, press, temp, sza, tauref, dist_sol,
     782
     783      subroutine phot(nlayer,
     784     $                lswitch, press, temp, sza, tauref, dist_sol,
    783785     $                rmco2, rmo3, j)
    784 c
     786
    785787c*****************************************************************
    786 c
     788
    787789      implicit none
    788 c
    789 #include "dimensions.h"
    790 #include "dimphys.h"
     790
     791!#include "dimensions.h"
     792!#include "dimphys.h"
    791793#include "chimiedata.h"
    792794#include "comcstfi.h"
    793 c
     795
    794796cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    795797c     inputs:
    796798cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    797 c       
    798       integer lswitch         ! interface level between chemistries
    799       real press(nlayermx)             ! pressure (hPa)
    800       real temp(nlayermx)              ! temperature (K)
    801       real sza                         ! solar zenith angle (deg)
    802       real tauref                      ! optical depth at 7 hpa
    803       real dist_sol                    ! sun distance (AU)
    804       real rmco2(nlayermx)             ! co2 mixing ratio
    805       real rmo3(nlayermx)              ! ozone mixing ratio
    806 c
     799       
     800      integer, intent(in) :: nlayer ! number of atmospheric layers
     801      integer :: lswitch            ! interface level between chemistries
     802      real :: press(nlayer)         ! pressure (hPa)
     803      real :: temp(nlayer)          ! temperature (K)
     804      real :: sza                   ! solar zenith angle (deg)
     805      real :: tauref                ! optical depth at 7 hpa
     806      real :: dist_sol              ! sun distance (AU)
     807      real :: rmco2(nlayer)         ! co2 mixing ratio
     808      real :: rmo3(nlayer)          ! ozone mixing ratio
     809
    807810cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    808811c     output:
    809812cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    810 c
    811       real j(nlayermx,nd)              ! interpolated photolysis rates (s-1)
    812 c
     813
     814      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1)
     815
    813816cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    814817c     local:
    815818cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    816 c
    817       integer icol, ij, indsza, indtau, indcol, indozo, indtemp,
    818      $        iozo, isza, itau, it, l
    819 c
    820       real col(nlayermx)               ! overhead air column   (molecule cm-2)
    821       real colo3(nlayermx)             ! overhead ozone column (molecule cm-2)
    822       real poids(2,2,2,2,2)            ! 5D interpolation weights
    823       real tref                        ! temperature  at 1.9 hPa in the gcm (K)
    824       real table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
    825       real cinf, csup, cicol,
    826      $     ciozo, cisza, citemp, citau
    827       real colo3min, dp, coef
    828       real ratio_o3(nlayermx)
    829       real tau
    830 c
     819
     820      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,
     821     $           iozo, isza, itau, it, l
     822
     823      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
     824      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
     825      real :: poids(2,2,2,2,2)            ! 5D interpolation weights
     826      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
     827      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
     828      real :: cinf, csup, cicol,
     829     $        ciozo, cisza, citemp, citau
     830      real :: colo3min, dp, coef
     831      real :: ratio_o3(nlayer)
     832      real :: tau
     833
    831834ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    832835c     day/night criterion
    833836ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    834 c
     837
    835838      if (sza .le. 95.) then
    836 c
     839
    837840ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    838841c     temperatures at 1.9 hPa in lookup table
    839842ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    840 c     
     843     
    841844      table_temp(1) = 226.2
    842845      table_temp(2) = 206.2
    843846      table_temp(3) = 186.2
    844847      table_temp(4) = 169.8
    845 c
     848
    846849ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    847850c     interpolation in solar zenith angle
    848851ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    849 c
     852
    850853      indsza = nsza - 1
    851854      do isza = 1,nsza
    … …  
    857860      cisza = (sza - szatab(indsza))
    858861     $       /(szatab(indsza + 1) - szatab(indsza))
    859 c
     862
    860863ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    861864c     interpolation in dust (tau)
    862865ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    863 c
     866
    864867      tau = min(tauref, tautab(ntau))
    865868      tau = max(tau, tautab(1))
    866 c
     869
    867870      indtau = ntau - 1
    868871      do itau = 1,ntau
    … …  
    874877      citau = (tau - tautab(indtau))
    875878     $       /(tautab(indtau + 1) - tautab(indtau))
    876 c
     879
    877880ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    878881c     co2 and ozone columns
    879882ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    880 c
     883
    881884c     co2 column at model top (molecule.cm-2)
    882 c
     885
    883886      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.
    884887     $                 /(mugaz*g)
    885 c
     888
    886889c     ozone column at model top
    887 c
     890
    888891      colo3(lswitch-1) = 0.
    889 c
     892
    890893c     co2 and ozone columns for other levels (molecule.cm-2)
    891 c
     894
    892895      do l = lswitch-2,1,-1
    893896         dp = (press(l) - press(l+1))*100.
    … …  
    900903     $              *6.022e22*dp/(mugaz*g)
    901904      end do
    902 c
     905
    903906c     ratio ozone column/minimal theoretical column (0.1 micron-atm)
    904907c
    … …  
    906909c     profile giving a column 0.1 micron-atmosphere at
    907910c     a surface pressure of 10 hpa.
    908 c
     911
    909912      do l = 1,lswitch-1
    910913         colo3min    = col(l)*7.171e-10
    … …  
    913916         ratio_o3(l) = max(ratio_o3(l), 1.)
    914917      end do
    915 c
     918
    916919ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    917920c     temperature dependence
    918921ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    919 c
     922
    920923c     1) search for temperature at 1.9 hPa (tref): vertical interpolation
    921 c
     924
    922925      tref = temp(1)
    923926      do l = (lswitch-1)-1,1,-1
    … …  
    931934      end do
    932935 10   continue
    933 c
     936
    934937c     2) interpolation in temperature
    935 c
     938
    936939      tref = min(tref,table_temp(1))
    937940      tref = max(tref,table_temp(ntemp))
    938 c
     941
    939942      do it = 2, ntemp
    940943         if (table_temp(it) .le. tref) then
    … …  
    946949      end do
    947950  20  continue
    948 c
     951
    949952ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    950953c     loop over vertical levels
    951954ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    952 c
     955
    953956      do l = 1,lswitch-1
    954 c
     957
    955958c     interpolation in air column
    956 c
     959
    957960         do icol = 0,200
    958961            if (colairtab(icol) .lt. col(l)) then
    … …  
    964967         end do
    965968 30      continue
    966 c
     969
    967970cc    interpolation in ozone column
    968 c
     971
    969972         indozo = nozo - 1
    970973         do iozo = 1,nozo
    … …  
    976979         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)
    977980     $          /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
    978 c
     981
    979982cc    4-dimensional interpolation weights
    980 c
     983
    981984c     poids(temp,sza,co2,o3,tau)
    982 c
     985
    983986         poids(1,1,1,1,1) = citemp
    984987     $                     *(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
    … …  
    10131016         poids(2,2,2,2,1) = (1.-citemp)
    10141017     $                     *cisza*(1.-cicol)*ciozo*(1.-citau)
    1015 c
     1018
    10161019         poids(1,1,1,1,2) = citemp
    10171020     $                     *(1.-cisza)*cicol*(1.-ciozo)*citau
    … …  
    10461049         poids(2,2,2,2,2) = (1.-citemp)
    10471050     $                     *cisza*(1.-cicol)*ciozo*citau
    1048 c
     1051
    10491052cc    4-dimensional interpolation in the lookup table
    1050 c
     1053
    10511054         do ij = 1,nd
    10521055           j(l,ij) =
    … …  
    10831086     $   + poids(2,2,2,2,1)
    10841087     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)
    1085 c
     1088
    10861089     $   + poids(1,1,1,1,2)
    10871090     $     *jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)
    … …  
    11171120     $     *jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
    11181121         end do
    1119 c
     1122
    11201123cc    correction for sun distance
    1121 c
     1124
    11221125         do ij = 1,nd
    11231126            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
    11241127         end do
    1125 c
     1128
    11261129ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11271130c     end of loop over vertical levels
    11281131ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1129 c
     1132
    11301133      end do
    1131 c
     1134
    11321135      else
    1133 c
     1136
    11341137ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11351138c     night
    11361139ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1137 c
     1140
    11381141      do ij = 1,nd
    11391142         do l = 1,lswitch-1
    … …  
    11411144         end do
    11421145      end do
    1143 c
     1146
    11441147      end if
    1145 c
     1148
    11461149      return
    11471150      end
    1148 c
     1151
    11491152c*****************************************************************
    1150 c
    1151       subroutine gcmtochim(zycol, lswitch, nesp, rm)
    1152 c
     1153
     1154      subroutine gcmtochim(nlayer, nq, zycol, lswitch, nesp, rm)
     1155
    11531156c*****************************************************************
    1154 c
    1155       use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d,
     1157
     1158      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
    11561159     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
    11571160     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
    11581161     &                      igcm_ch4
     1162
    11591163      implicit none
    1160 c
    1161 #include "dimensions.h"
    1162 #include "dimphys.h"
     1164
     1165!#include "dimensions.h"
     1166!#include "dimphys.h"
    11631167#include "callkeys.h"
    11641168!#include "tracer.h"
    1165 c
     1169
    11661170cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11671171c     inputs:
    11681172cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1169 c     
    1170       real zycol(nlayermx,nqmx)! species volume mixing ratio in the gcm
    1171 c
    1172       integer nesp             ! number of species in the chemistry
    1173       integer lswitch          ! interface level between chemistries
    1174 c
     1173     
     1174      integer, intent(in) :: nlayer ! number of atmospheric layers
     1175      integer, intent(in) :: nq     ! number of tracers
     1176      real :: zycol(nlayer,nq)      ! species volume mixing ratio in the gcm
     1177
     1178      integer :: nesp               ! number of species in the chemistry
     1179      integer :: lswitch            ! interface level between chemistries
     1180
    11751181cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11761182c     outputs:
    11771183cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1178 c
    1179       real rm(nlayermx,nesp)   ! species volume mixing ratio
    1180 c     
     1184
     1185      real :: rm(nlayer,nesp)   ! species volume mixing ratio
     1186     
    11811187cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    11821188c     local:
    11831189cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1184 c
    1185       integer l,iq
     1190
     1191      integer :: l, iq
    11861192     
    11871193c     tracer indexes in the chemistry:
    … …  
    12031209      integer,parameter :: i_hox  = 15
    12041210      integer,parameter :: i_ox   = 16
    1205 c     
     1211
    12061212cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12071213c     initialise chemical species
    12081214cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1209 c
     1215
    12101216      do l = 1,lswitch-1
    12111217         rm(l,i_co2)  = max(zycol(l, igcm_co2),     1.e-30)
    … …  
    12331239         end do
    12341240      end if
    1235 c
     1241
    12361242cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12371243c     initialise chemical families                     c
    12381244cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1239 c
     1245
    12401246      do l = 1,lswitch-1
    12411247         rm(l,i_hox) = rm(l,i_h)
    … …  
    12451251     $               + rm(l,i_o3)
    12461252      end do
    1247 c
     1253
    12481254      return
    12491255      end
    1250 c
     1256
    12511257c*****************************************************************
    12521258c
    1253       subroutine chimtogcm(zycol, lswitch, nesp, rm)
     1259      subroutine chimtogcm(nlayer, nq, zycol, lswitch, nesp, rm)
    12541260c
    12551261c*****************************************************************
    12561262c
    1257       use tracer_mod, only: nqmx, igcm_co2, igcm_co, igcm_o, igcm_o1d,
     1263      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
    12581264     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
    12591265     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_h2o_vap,
    12601266     &                      igcm_ch4
     1267
    12611268      implicit none
    1262 c
    1263 #include "dimensions.h"
    1264 #include "dimphys.h"
     1269
     1270!#include "dimensions.h"
     1271!#include "dimphys.h"
    12651272#include "callkeys.h"
    12661273!#include "tracer.h"
    1267 c
     1274
    12681275cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12691276c     inputs:
    12701277cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1271 c
    1272       integer nesp               ! number of species in the chemistry
    1273       integer lswitch            ! interface level between chemistries
    1274 c     
    1275       real rm(nlayermx,nesp)     ! species volume mixing ratio
    1276 c
     1278
     1279      integer, intent(in) :: nlayer ! number of atmospheric layers
     1280      integer, intent(in) :: nq     ! number of tracers
     1281      integer :: nesp               ! number of species in the chemistry
     1282      integer :: lswitch            ! interface level between chemistries
     1283     
     1284      real :: rm(nlayer,nesp)       ! species volume mixing ratio
     1285
    12771286cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12781287c     output:
    12791288cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1280 c     
    1281       real zycol(nlayermx,nqmx)  ! species volume mixing ratio in the gcm
    1282 c
     1289
     1290      real :: zycol(nlayer,nq)  ! species volume mixing ratio in the gcm
     1291
    12831292cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    12841293c     local:
    12851294cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1286 c
    1287       integer l,iq
     1295
     1296      integer l, iq
    12881297     
    12891298c     tracer indexes in the chemistry:
    … …  
    13051314      integer,parameter :: i_hox  = 15
    13061315      integer,parameter :: i_ox   = 16
    1307 c     
     1316
    13081317cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    13091318c     save mixing ratios for the gcm
    13101319cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1311 c
     1320
    13121321      do l = 1,lswitch-1
    13131322         zycol(l, igcm_co2)     = rm(l,i_co2)
    … …  
    13311340         end do
    13321341      end if
    1333 c
     1342
    13341343      return
    13351344      end
    1336 c
     1345
    13371346c*****************************************************************
    1338 c
    1339       subroutine chemrates(lswitch, dens, press, t,
     1347
     1348      subroutine chemrates(nlayer,
     1349     $                     lswitch, dens, press, t,
    13401350     $                     surfdust1d, surfice1d,
    13411351     $                     a001, a002, a003,
    … …  
    13491359     $                     h001, h002, h003, h004, h005,
    13501360     $                     t001, t002, t003, tau)
    1351 c
     1361
    13521362c*****************************************************************
    1353 c
     1363
    13541364      implicit none
    1355 c
    1356 #include "dimensions.h"
    1357 #include "dimphys.h"
     1365
     1366!#include "dimensions.h"
     1367!#include "dimphys.h"
    13581368#include "comcstfi.h"
    1359 c
     1369
    13601370cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    13611371c     inputs:                                                    c
    13621372cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1363 c
    1364       integer lswitch                  ! interface level between chemistries
    1365 
    1366       real dens(nlayermx)              ! density (cm-3)
    1367       real press(nlayermx)             ! pressure (hpa)
    1368       real t(nlayermx)                 ! temperature (k)
    1369       real surfdust1d(nlayermx)        ! dust surface area (cm^2/cm^3)
    1370       real surfice1d(nlayermx)         ! ice surface area (cm^2/cm^3)
    1371       real tribo                       ! switch for triboelectricity
    1372       real tau                         ! dust opacity at 7 hpa
    1373 c
    1374       parameter (tribo   = 0.)        ! switch for triboelectricity
    1375 c
     1373
     1374      integer, intent(in) :: nlayer   ! number of atmospheric layers
     1375      integer :: lswitch              ! interface level between chemistries
     1376
     1377      real :: dens(nlayer)            ! density (cm-3)
     1378      real :: press(nlayer)           ! pressure (hpa)
     1379      real :: t(nlayer)               ! temperature (k)
     1380      real :: surfdust1d(nlayer)      ! dust surface area (cm^2/cm^3)
     1381      real :: surfice1d(nlayer)       ! ice surface area (cm^2/cm^3)
     1382      real :: tau                     ! dust opacity at 7 hpa
     1383
     1384      real, parameter :: tribo   = 0. ! switch for triboelectricity
     1385
    13761386cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    13771387c     outputs:                                                   c
    13781388cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1379 c
    1380       real a001(nlayermx), a002(nlayermx), a003(nlayermx)
    1381       real b001(nlayermx), b002(nlayermx), b003(nlayermx),
    1382      $     b004(nlayermx), b005(nlayermx), b006(nlayermx),
    1383      $     b007(nlayermx), b008(nlayermx), b009(nlayermx)
    1384       real c001(nlayermx), c002(nlayermx), c003(nlayermx),
    1385      $     c004(nlayermx), c005(nlayermx), c006(nlayermx),
    1386      $     c007(nlayermx), c008(nlayermx), c009(nlayermx),
    1387      $     c010(nlayermx), c011(nlayermx), c012(nlayermx),
    1388      $     c013(nlayermx), c014(nlayermx), c015(nlayermx),
    1389      $     c016(nlayermx), c017(nlayermx), c018(nlayermx)
    1390       real d001(nlayermx), d002(nlayermx), d003(nlayermx)
    1391       real e001(nlayermx), e002(nlayermx), e003(nlayermx)
    1392       real h001(nlayermx), h002(nlayermx), h003(nlayermx),
    1393      $     h004(nlayermx), h005(nlayermx)
    1394       real t001(nlayermx), t002(nlayermx), t003(nlayermx)
    1395 c
     1389
     1390      real :: a001(nlayer), a002(nlayer), a003(nlayer)
     1391      real :: b001(nlayer), b002(nlayer), b003(nlayer),
     1392     $        b004(nlayer), b005(nlayer), b006(nlayer),
     1393     $        b007(nlayer), b008(nlayer), b009(nlayer)
     1394      real :: c001(nlayer), c002(nlayer), c003(nlayer),
     1395     $        c004(nlayer), c005(nlayer), c006(nlayer),
     1396     $        c007(nlayer), c008(nlayer), c009(nlayer),
     1397     $        c010(nlayer), c011(nlayer), c012(nlayer),
     1398     $        c013(nlayer), c014(nlayer), c015(nlayer),
     1399     $        c016(nlayer), c017(nlayer), c018(nlayer)
     1400      real :: d001(nlayer), d002(nlayer), d003(nlayer)
     1401      real :: e001(nlayer), e002(nlayer), e003(nlayer)
     1402      real :: h001(nlayer), h002(nlayer), h003(nlayer),
     1403     $        h004(nlayer), h005(nlayer)
     1404      real :: t001(nlayer), t002(nlayer), t003(nlayer)
     1405
    13961406cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    13971407c     local:                                                     c
    13981408cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1399 c
    1400       real ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
    1401       real ef, efmax, lossh2o, lossch4, lossco2
    1402 c
    1403       integer l
    1404       real k1a, k1b, k1a0, k1b0, k1ainf
    1405       real x, y, fc, fx
    1406 c
     1409
     1410      real :: ak0, ak1, rate, rate1, rate2, xpo, xpo1, xpo2
     1411      real :: ef, efmax, lossh2o, lossch4, lossco2
     1412
     1413      integer :: l
     1414      real :: k1a, k1b, k1a0, k1b0, k1ainf
     1415      real :: x, y, fc, fx
     1416
    14071417cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    14081418c     compute reaction rates
    14091419cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1410 c
     1420
    14111421      do l = 1,lswitch-1
    1412 c
     1422
    14131423cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    14141424c        oxygen compounds
Note: See TracChangeset for help on using the changeset viewer.