source: trunk/LMDZ.MARS/libf/aeronomars/calchim.F @ 551

      subroutine calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
     $                   zzlev,zzlay,zday,pq,pdq,dqchim,dqschim,dqcloud,
     $                   dqscloud,tauref,co2ice,
     $                   pu,pdu,pv,pdv,surfdust,surfice)
c
      implicit none
c
c=======================================================================
c
c   subject:
c   --------
c
c  Prepare the call for the photochemical module, and send back the
c  tendencies from photochemistry in the chemical species mass mixing ratios
c
c   Author:   Sebastien Lebonnois (08/11/2002)
c   -------
c    update 12/06/2003 for water ice clouds and compatibility with dust
c    update 07/2003 for coupling with thermosphere (Monica Angelats-i-Coll)
c    update 03/05/2005 cosmetic changes (Franck Lefevre)
c    update sept. 2008 identify tracers by their names (Ehouarn Millour)
c    update 17/03/2011 synchronize with latest version of chemistry (Franck Lefevre)
c    update 05/12/2011 synchronize with latest version of chemistry (Franck Lefevre)
c
c   Arguments:
c   ----------
c
c  Input:
c
c    ptimestep                  timestep (s)
c    pplay(ngridmx,nlayermx)    Pressure at the middle of the layers (Pa)
c    pplev(ngridmx,nlayermx+1)  Intermediate pressure levels (Pa)
c    pt(ngridmx,nlayermx)       Temperature (K)
c    pdt(ngridmx,nlayermx)      Temperature tendency (K)
c    pu(ngridmx,nlayermx)       u component of the wind (ms-1)
c    pdu(ngridmx,nlayermx)      u component tendency (K)
c    pv(ngridmx,nlayermx)       v component of the wind (ms-1)
c    pdv(ngridmx,nlayermx)      v component tendency (K)
c    dist_sol                   distance of the sun (AU)
c    mu0(ngridmx)               cos of solar zenith angle (=1 when sun at zenith)
c    pq(ngridmx,nlayermx,nqmx)  Advected fields, ie chemical species here
c    pdq(ngridmx,nlayermx,nqmx) Previous tendencies on pq
c    tauref(ngridmx)            Optical depth at 7 hPa
c    co2ice(ngridmx)            co2 ice surface layer (kg.m-2)
c    surfdust(ngridmx,nlayermx) dust surface area (m2/m3)
c    surfice(ngridmx,nlayermx)  ice surface area (m2/m3)
c
c  Output:
c
c    dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
c    dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
c
c=======================================================================

c    Declarations :
c    --------------

#include "dimensions.h"
#include "dimphys.h"
#include "chimiedata.h"
#include "tracer.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "conc.h"

c Arguments :
c -----------

c   inputs:
c   -------

      real    ptimestep
      real    pplay(ngridmx,nlayermx)    ! pressure at the middle of the layers
      real    zzlay(ngridmx,nlayermx)    ! pressure at the middle of the layers
      real    pplev(ngridmx,nlayermx+1)  ! intermediate pressure levels
      real    zzlev(ngridmx,nlayermx+1)  ! altitude at layer boundaries
      real    pt(ngridmx,nlayermx)       ! temperature
      real    pdt(ngridmx,nlayermx)      ! temperature tendency
      real    pu(ngridmx,nlayermx)       ! u component of the wind (m.s-1)
      real    pdu(ngridmx,nlayermx)      ! u component tendency
      real    pv(ngridmx,nlayermx)       ! v component of the wind (m.s-1)
      real    pdv(ngridmx,nlayermx)      ! v component tendency
      real    dist_sol                   ! distance of the sun (AU)
      real    mu0(ngridmx)               ! cos of solar zenith angle (=1 when sun at zenith)
      real    pq(ngridmx,nlayermx,nqmx)  ! tracers mass mixing ratio
      real    pdq(ngridmx,nlayermx,nqmx) ! previous tendencies
      real    zday                       ! date (time since Ls=0, in martian days)
      real    tauref(ngridmx)            ! optical depth at 7 hPa
      real    co2ice(ngridmx)            ! co2 ice surface layer (kg.m-2)
      real    surfdust(ngridmx,nlayermx) ! dust surface area (m2/m3)
      real    surfice(ngridmx,nlayermx)  !  ice surface area (m2/m3)

c   outputs:
c   --------

      real dqchim(ngridmx,nlayermx,nqmx) ! tendencies on pq due to chemistry
      real dqschim(ngridmx,nqmx)         ! tendencies on qsurf 
      real dqcloud(ngridmx,nlayermx,nqmx)! tendencies on pq due to condensation
      real dqscloud(ngridmx,nqmx)        ! tendencies on qsurf 

c Local variables :
c -----------------

      character*20 str20
      integer  ig,l,i,iq
      integer  foundswitch, lswitch
      integer  ig_vl1
      real latvl1, lonvl1
      real zq(ngridmx,nlayermx,nqmx) ! pq+pdq*ptimestep before chemistry
                                     ! new mole fraction after
      real zt(ngridmx,nlayermx)      ! temperature
      real zu(ngridmx,nlayermx)      ! u component of the wind
      real zv(ngridmx,nlayermx)      ! v component of the wind
      real taucol                    ! optical depth at 7 hPa

      logical depos                  ! switch for dry deposition
      parameter (depos=.false.)
c
c  for each column of atmosphere:
c
      real zpress(nlayermx)       !  Pressure (mbar)
      real zdens(nlayermx)        !  Density  (cm-3)
      real ztemp(nlayermx)        !  Temperature (K)
      real zlocal(nlayermx)       !  Altitude (km)
      real zycol(nlayermx,nqmx)   !  Composition (mole fractions)
      real szacol                 !  Solar zenith angle
      real surfice1d(nlayermx)    !  Ice surface area (cm2/cm3)
      real surfdust1d(nlayermx)   !  Dust surface area (cm2/cm3)
      real jo3(nlayermx)          !  Photodissociation rate O3->O1D (s-1)
c
c for output:
c
      real jo3_3d(ngridmx,nlayermx)  !  Photodissociation rate O3->O1D (s-1)

      logical output              ! to issue calls to writediagfi and stats
      parameter (output=.true.)   ! see at end of routine

      logical,save :: firstcall=.true.
      integer,save :: nbq       ! number of tracers used in the chemistry
      integer,save :: niq(nqmx) ! array storing the indexes of the tracers

! index of tracers:

      integer,save :: i_co2=0
      integer,save :: i_co=0
      integer,save :: i_o=0
      integer,save :: i_o1d=0
      integer,save :: i_o2=0
      integer,save :: i_o3=0
      integer,save :: i_h=0
      integer,save :: i_h2=0
      integer,save :: i_oh=0
      integer,save :: i_ho2=0
      integer,save :: i_h2o2=0
      integer,save :: i_ch4=0
      integer,save :: i_n2=0
      integer,save :: i_ar=0
      integer,save :: i_ice=0 ! water ice
      integer,save :: i_h2o=0 ! water vapour
c
c=======================================================================
c     initialization of the chemistry (first call only)
c=======================================================================
c
      if (firstcall) then
c
         if (photochem) then
            print*,'calchim: Read photolysis lookup table'
            call read_phototable
         end if

         ! find index of chemical tracers to use
         nbq=0 ! to count number of tracers
         i_co2=igcm_co2
         if (i_co2.eq.0) then
           write(*,*) "calchim: Error; no CO2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_co2
         endif
         i_co=igcm_co
         if (i_co.eq.0) then
           write(*,*) "calchim: Error; no CO tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_co
         endif
         i_o=igcm_o
         if (i_o.eq.0) then
           write(*,*) "calchim: Error; no O tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o
         endif
         i_o1d=igcm_o1d
         if (i_o1d.eq.0) then
           write(*,*) "calchim: Error; no O1D tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o1d
         endif
         i_o2=igcm_o2
         if (i_o2.eq.0) then
           write(*,*) "calchim: Error; no O2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o2
         endif
         i_o3=igcm_o3
         if (i_o3.eq.0) then
           write(*,*) "calchim: Error; no O3 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_o3
         endif
         i_h=igcm_h
         if (i_h.eq.0) then
           write(*,*) "calchim: Error; no H tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h
         endif
         i_h2=igcm_h2
         if (i_h2.eq.0) then
           write(*,*) "calchim: Error; no H2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2
         endif
         i_oh=igcm_oh
         if (i_oh.eq.0) then
           write(*,*) "calchim: Error; no OH tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_oh
         endif
         i_ho2=igcm_ho2
         if (i_ho2.eq.0) then
           write(*,*) "calchim: Error; no HO2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ho2
         endif
         i_h2o2=igcm_h2o2
         if (i_h2o2.eq.0) then
           write(*,*) "calchim: Error; no H2O2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2o2
         endif
         i_ch4=igcm_ch4
         if (i_ch4.eq.0) then
           write(*,*) "calchim: Error; no CH4 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ch4
         endif
         i_n2=igcm_n2
         if (i_n2.eq.0) then
           write(*,*) "calchim: Error; no N2 tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_n2
         endif
         i_ar=igcm_ar
         if (i_ar.eq.0) then
           write(*,*) "calchim: Error; no AR tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ar
         endif
         i_ice=igcm_h2o_ice
         if (i_ice.eq.0) then
           write(*,*) "calchim: Error; no water ice tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_ice
         endif
         i_h2o=igcm_h2o_vap
         if (i_h2o.eq.0) then
           write(*,*) "calchim: Error; no water vapor tracer !!!"
           stop
         else
           nbq=nbq+1
           niq(nbq)=i_h2o
         endif

         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
         write(*,*) '               i_co2  = ',i_co2
         write(*,*) '               i_co   = ',i_co
         write(*,*) '               i_o    = ',i_o
         write(*,*) '               i_o1d  = ',i_o1d
         write(*,*) '               i_o2   = ',i_o2
         write(*,*) '               i_o3   = ',i_o3
         write(*,*) '               i_h    = ',i_h
         write(*,*) '               i_h2   = ',i_h2
         write(*,*) '               i_oh   = ',i_oh
         write(*,*) '               i_ho2  = ',i_ho2
         write(*,*) '               i_h2o2 = ',i_h2o2
         write(*,*) '               i_ch4  = ',i_ch4
         write(*,*) '               i_n2   = ',i_n2
         write(*,*) '               i_ar   = ',i_ar
         write(*,*) '               i_ice  = ',i_ice
         write(*,*) '               i_h2o  = ',i_h2o
         
         firstcall = .false.
      end if ! if (firstcall)

! Initialize output tendencies to zero (to handle case of tracers which
! are not used in the chemistry (e.g. dust))

      dqchim(:,:,:) = 0
      dqschim(:,:)  = 0
c
!     latvl1= 22.27
!     lonvl1= -47.94
!     ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.)  -2 )*iim +
!    $        int(1.5+(lonvl1+180)*iim/360.)
c
c=======================================================================
c     loop over grid
c=======================================================================
c
      do ig = 1,ngridmx
c
c     local updates
c 
         foundswitch = 0
         do l = 1,nlayermx
            do i = 1,nbq
              iq = niq(i) ! get tracer index
              zq(ig,l,iq) = pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep
              zycol(l,iq) = zq(ig,l,iq)*mmean(ig,l)/mmol(iq)
            end do
            zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
            zu(ig,l) = pu(ig,l) + pdu(ig,l)*ptimestep
            zv(ig,l) = pv(ig,l) + pdv(ig,l)*ptimestep
c
            zpress(l) = pplay(ig,l)/100.
            ztemp(l)  = zt(ig,l)
            zdens(l)  = zpress(l)/(kb*1.e4*ztemp(l))
            zlocal(l) = zzlay(ig,l)/1000.
c
c           surfdust1d and surfice1d: conversion from m2/m3 to cm2/cm3
c
            surfdust1d(l) = surfdust(ig,l)*1.e-2
            surfice1d(l)  = surfice(ig,l)*1.e-2
c
c           search for switch index between regions
c
            if (photochem .and. thermochem) then
               if (foundswitch .eq. 0 .and. pplay(ig,l).lt.1.e-3) then
                  lswitch = l
                  foundswitch=1
               end if
            end if
            if ( .not. photochem) then
               lswitch = 22
            end if
            if (.not.  thermochem) then
               lswitch = min(50,nlayermx+1)    ! FL (original value: 33)
            end if
c
         end do ! of do l=1,nlayermx
c
         szacol = acos(mu0(ig))*180./pi
         taucol = tauref(ig)
c
c=======================================================================
c     call chemical subroutine
c=======================================================================
c
c        chemistry in lower atmosphere
c
         if (photochem) then
            call photochemistry(lswitch,zycol,szacol,ptimestep,
     $                          zpress,ztemp,zdens,dist_sol,
     $                          surfdust1d,surfice1d,jo3,taucol)
         end if
c
c        chemistry in upper atmosphere
c
         if (thermochem) then
            call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,
     $                       zlocal,szacol,ptimestep,zday)
         end if
c
c        dry deposition
c
         if (depos) then
            call deposition(ig, ig_vl1, pplay, pplev, zzlay, zzlev, 
     $                      zu, zv, zt, zycol, ptimestep, co2ice)
         end if
c
c=======================================================================
c     tendencies
c=======================================================================
c
c     must be 0. for water ice:
c
         if (water) then
            do l = 1,nlayermx
               dqchim(ig,l,i_ice) = 0.
            end do
         end if
c
c     tendency for CO2 = - sum of others for lower atmosphere
c     tendency for O   = - sum of others for upper atmosphere
c
         do l = 1,nlayermx
            if (l .lt. lswitch) then
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                  if ((iq.ne.i_co2).and.(iq.ne.i_ice)) then
                     dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)/mmean(ig,l)
     $                                - zq(ig,l,iq))/ptimestep
                  else if (iq.eq.i_co2) then
                     dqchim(ig,l,iq) = 0.
                  end if
                  dqschim(ig,iq) = 0.
               end do ! of do i=1,nbq
               
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                  if (iq.ne.i_co2) then
                     dqchim(ig,l,i_co2) = dqchim(ig,l,i_co2) 
     $                                  - dqchim(ig,l,iq)
                  end if
               end do
             else if (l .ge. lswitch) then
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                   if ((iq.ne.i_o).and.(iq.ne.i_ice)) then
                      dqchim(ig,l,iq) = (zycol(l,iq)*mmol(iq)
     $                                  /mmean(ig,l)
     $                                 - zq(ig,l,iq))/ptimestep
                   else if (iq.eq.i_o) then
                      dqchim(ig,l,iq) = 0.
                   end if
               enddo
               
               do i=1,nbq
                 iq=niq(i) ! get tracer index
                   if (iq.ne.i_o) then
                      dqchim(ig,l,i_o) = dqchim(ig,l,i_o) 
     $                                 - dqchim(ig,l,iq)
                   end if
               end do
             end if ! of if (l.lt.lswitch) else if (l.ge.lswitch)
          end do ! of do l = 1,nlayermx
c
c     condensation of h2o2
c
          call perosat(ig,ptimestep,pplev,pplay,
     $                 ztemp,zycol,dqcloud,dqscloud)
c
c     density and j(o3->o1d), for outputs
c
          do l = 1,nlayermx
             jo3_3d(ig,l) = jo3(l)
          end do
c
c=======================================================================
c     end of loop over grid
c=======================================================================
c
      end do ! of do ig=1,ngridmx
c
c=======================================================================
c     write outputs
c=======================================================================
c
! value of parameter 'output' to trigger writting of outputs
! is set above at the declaration of the variable.

      if (output) then
         if (ngridmx .gt. 1) then
            call writediagfi(ngridmx,'jo3','j o3->o1d',
     $                       's-1',3,jo3_3d(1,1))
           if (callstats) then
            call wstats(ngridmx,'jo3','j o3->o1d',
     $                       's-1',3,jo3_3d(1,1))
           endif
         end if ! of if (ngridmx.gt.1)
      endif ! of if (output)
c
      end