source: trunk/LMDZ.MARS/libf/aeronomars/inichim_newstart.F90 @ 635

      subroutine inichim_newstart(pq, qsurf, ps, flagh2o, flagthermo)

      implicit none

!=======================================================================
!
!  Purpose:
!  --------
!
!  Initialization of the chemistry for use in the building of a new start file
!  used by program newstart.F
!  also used by program testphys1d.F
!
!  Authors: 
!  -------
!  Initial version 11/2002 by Sebastien Lebonnois
!  Revised 07/2003 by Monica Angelats-i-Coll to use input files
!  Modified 10/2008 Identify tracers by their names Ehouarn Millour
!  Modified 11/2011 Addition of methane Franck Lefevre
!  Rewritten 04/2012 Franck Lefevre
!
!  Arguments:
!  ----------
!
!  pq(iip1,jjp1,llm,nqmx)  Advected fields, ie chemical species here
!  qsurf(ngridmx,nqmx)     Amount of tracer on the surface (kg/m2)
!  ps(iip1,jjp1)           Surface pressure (Pa)
!  flagh2o                 flag for initialisation of h2o (1: yes / 0: no)
!  flagthermo              flag for initialisation of thermosphere only (1: yes / 0: no)
!
!=======================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "paramet.h"
#include "tracer.h"
#include "comvert.h"
#include "callkeys.h"
#include "datafile.h"

! inputs :

      real,intent(in) :: ps(iip1,jjp1)            ! surface pressure in the gcm (Pa)   
      integer,intent(in) :: flagh2o               ! flag for h2o initialisation
      integer,intent(in) :: flagthermo            ! flag for thermosphere initialisation only

! outputs :

      real,intent(out) :: pq(iip1,jjp1,llm,nqmx)  ! advected fields, ie chemical species
      real,intent(out) :: qsurf(ngridmx,nqmx)     ! surface values (kg/m2) of tracers

! local :

      integer :: iq, i, j, l, n, nbqchem
      integer :: count, ierr, dummy
      real    :: mmean(iip1,jjp1,llm)             ! mean molecular mass (g)
      real    :: pgcm                             ! pressure at each layer in the gcm (Pa)

      integer, parameter         :: nalt = 252    ! number of altitudes in the initialization files
      integer, parameter         :: nspe = 14     ! number of species in the initialization files
      integer, dimension(nspe)   :: niq           ! local index of species in initialization files
      real, dimension(nalt)      :: tinit, zzfile ! temperature in initialization files
      real, dimension(nalt)      :: pinit         ! pressure in initialization files
      real, dimension(nalt)      :: densinit      ! total number density in initialization files
      real, dimension(nalt,nspe) :: vmrinit       ! mixing ratios in initialization files
      real, dimension(nspe)      :: vmrint        ! mixing ratio interpolated onto the gcm vertical grid
      real                       :: vmr

      character(len=20)          :: txt           ! to store some text

! 1. identify tracers by their names: (and set corresponding values of mmol)

! 1.1 initialize tracer indexes to zero:

      do iq = 1,nqmx
        igcm_dustbin(iq) = 0
      end do

      igcm_dust_mass   = 0
      igcm_dust_number = 0
      igcm_ccn_mass    = 0
      igcm_ccn_number  = 0
      igcm_h2o_vap     = 0
      igcm_h2o_ice     = 0
      igcm_co2         = 0
      igcm_co          = 0
      igcm_o           = 0
      igcm_o1d         = 0
      igcm_o2          = 0
      igcm_o3          = 0
      igcm_h           = 0
      igcm_h2          = 0
      igcm_oh          = 0
      igcm_ho2         = 0
      igcm_h2o2        = 0
      igcm_ch4         = 0
      igcm_n2          = 0
      igcm_ar          = 0
      igcm_ar_n2       = 0
      igcm_co2plus     = 0
      igcm_oplus       = 0
      igcm_o2plus      = 0
      igcm_coplus      = 0
      igcm_cplus       = 0
      igcm_nplus       = 0
      igcm_noplus      = 0
      igcm_n2plus      = 0
      igcm_hplus       = 0
      igcm_elec        = 0

! 1.2 find dust tracers

      count = 0

      if (dustbin > 0) then
         do iq = 1,nqmx
            txt = " "
            write(txt,'(a4,i2.2)') 'dust', count + 1
            if (noms(iq) == txt) then
               count = count + 1
               igcm_dustbin(count) = iq
               mmol(iq) = 100.
            end if
         end do !do iq=1,nqmx
      end if ! of if (dustbin.gt.0)

      if (doubleq) then
         do iq = 1,nqmx
            if (noms(iq) == "dust_mass") then
               igcm_dust_mass = iq
               count = count + 1
            end if
            if (noms(iq) == "dust_number") then
               igcm_dust_number = iq
               count = count + 1
            end if
         end do
      end if ! of if (doubleq)

      if (scavenging) then
         do iq = 1,nqmx
            if (noms(iq) == "ccn_mass") then
               igcm_ccn_mass = iq
               count = count + 1
            end if
            if (noms(iq) == "ccn_number") then
               igcm_ccn_number = iq
               count = count + 1
            end if
         end do
      end if ! of if (scavenging)

      if (submicron) then
         do iq=1,nqmx
            if (noms(iq) == "dust_submicron") then
               igcm_dust_submicron = iq
               mmol(iq) = 100.
               count = count + 1
            end if
         end do
      end if ! of if (submicron)

! 1.3 find chemistry and water tracers

      nbqchem = 0
      do iq = 1,nqmx
         if (noms(iq) == "co2") then
            igcm_co2 = iq
            mmol(igcm_co2) = 44.
            count = count + 1
            nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "co") then
           igcm_co = iq
           mmol(igcm_co) = 28.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "o") then
           igcm_o = iq
           mmol(igcm_o) = 16.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "o1d") then
           igcm_o1d = iq
           mmol(igcm_o1d) = 16.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "o2") then
           igcm_o2 = iq
           mmol(igcm_o2) = 32.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "o3") then
           igcm_o3 = iq
           mmol(igcm_o3) = 48.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "h") then
           igcm_h = iq
           mmol(igcm_h) = 1.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "h2") then
           igcm_h2 = iq
           mmol(igcm_h2) = 2.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "oh") then
           igcm_oh = iq
           mmol(igcm_oh) = 17.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "ho2") then
           igcm_ho2 = iq
           mmol(igcm_ho2) = 33.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "h2o2") then
           igcm_h2o2 = iq
           mmol(igcm_h2o2) = 34.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "ch4") then
           igcm_ch4 = iq
           mmol(igcm_ch4) = 16.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "n2") then
           igcm_n2 = iq
           mmol(igcm_n2) = 28.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "n") then
           igcm_n = iq
           mmol(igcm_n) = 14.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "n2d") then
           igcm_n2d = iq
           mmol(igcm_n2d) = 14.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "no") then
           igcm_no = iq
           mmol(igcm_no) = 30.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "no2") then
           igcm_no2 = iq
           mmol(igcm_no2) = 46.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "ar") then
           igcm_ar = iq
           mmol(igcm_ar) = 40.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "h2o_vap") then
           igcm_h2o_vap = iq
           mmol(igcm_h2o_vap) = 18.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "h2o_ice") then
           igcm_h2o_ice = iq
           mmol(igcm_h2o_ice) = 18.
           count = count + 1
           nbqchem = nbqchem + 1
        end if

! 1.4 find ions

        if (noms(iq) == "co2plus") then
           igcm_co2plus = iq
           mmol(igcm_co2plus) = 44.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "oplus") then
           igcm_oplus = iq
           mmol(igcm_oplus) = 16.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "o2plus") then
           igcm_o2plus = iq
           mmol(igcm_o2plus) = 32.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "coplus") then
           igcm_coplus = iq
           mmol(igcm_coplus) = 28.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "cplus") then
           igcm_cplus = iq
           mmol(igcm_cplus) = 12.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "nplus") then
           igcm_nplus = iq
           mmol(igcm_nplus) = 14.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "noplus") then
           igcm_noplus = iq
           mmol(igcm_noplus) = 30.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "n2plus") then
           igcm_n2plus = iq
           mmol(igcm_n2plus) = 28.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "hplus") then
           igcm_hplus = iq
           mmol(igcm_hplus) = 1.
           count = count + 1
           nbqchem = nbqchem + 1
        end if
        if (noms(iq) == "elec") then
           igcm_elec = iq
           mmol(igcm_elec) = 1./1822.89
           count = count + 1
        end if

! 1.5 find idealized non-condensible tracer

        if (noms(iq) == "Ar_N2") then
           igcm_ar_n2 = iq
           mmol(igcm_ar_n2) = 30.
           count = count + 1
        end if

      end do ! of do iq=1,nqmx
      
! 1.6 check that we identified all tracers:

      if (count /= nqmx) then
         write(*,*) "inichim_newstart: found only ",count," tracers"
         write(*,*) "                  expected ",nqmx
         do iq = 1,count
            write(*,*) '      ', iq, ' ', trim(noms(iq))
         end do
         stop
      else
         write(*,*) "inichim_newstart: found all expected tracers"
         do iq = 1,nqmx
            write(*,*) '      ', iq, ' ', trim(noms(iq))
         end do
      end if

! 2. load in chemistry data for initialization:

! order of major species in initialization file:
!
!    1: co2 
!    2: ar
!    3: n2  
!    4: o2  
!    5: co  
!    6: o   
!    7: h2
!
! order of minor species in initialization file:
!
!    1: h  
!    2: oh 
!    3: ho2 
!    4: h2o 
!    5: h2o2 
!    6: o1d
!    7: o3

! major species:

      niq(1) = igcm_co2
      niq(2) = igcm_ar
      niq(3) = igcm_n2
      niq(4) = igcm_o2
      niq(5) = igcm_co
      niq(6) = igcm_o
      niq(7) = igcm_h2

! minor species:

      niq(8)  = igcm_h
      niq(9)  = igcm_oh
      niq(10) = igcm_ho2
      niq(11) = igcm_h2o_vap
      niq(12) = igcm_h2o2
      niq(13) = igcm_o1d 
      niq(14) = igcm_o3

! 2.1 open initialization files

      open(210, iostat=ierr,file=trim(datafile)// &
                            '/atmosfera_LMD_may.dat')
      if (ierr /= 0) then
         write(*,*)'Error : cannot open file atmosfera_LMD_may.dat '
         write(*,*)'(in aeronomars/inichim.F)'
         write(*,*)'It should be in :', trim(datafile),'/'
         write(*,*)'1) You can change this path in callphys.def with'
         write(*,*)'   datadir=/path/to/datafiles/'
         write(*,*)'2) If necessary atmosfera_LMD_may.dat (and others)'
         write(*,*)'   can be obtained online on:'
         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
         stop
      end if
      open(220, iostat=ierr,file=trim(datafile)// &
                            '/atmosfera_LMD_min.dat')
      if (ierr /= 0) then
         write(*,*)'Error : cannot open file atmosfera_LMD_min.dat '
         write(*,*)'(in aeronomars/inichim.F)'
         write(*,*)'It should be in :', trim(datafile),'/'
         write(*,*)'1) You can change this path in callphys.def with'
         write(*,*)'   datadir=/path/to/datafiles/'
         write(*,*)'2) If necessary atmosfera_LMD_min.dat (and others)'
         write(*,*)'   can be obtained online on:'
         write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
         stop
      end if

! 2.2 read initialization files

! major species

      read(210,*)
      do l = 1,nalt
         read(210,*) dummy, tinit(l), pinit(l), densinit(l), &
                     (vmrinit(l,n), n = 1,7)
         pinit(l) = pinit(l)*100.              ! conversion in Pa
         pinit(l) = log(pinit(l))              ! for the vertical interpolation
      end do
      close(210)

! minor species

      read(220,*)
      do l = 1,nalt
         read(220,*) dummy, (vmrinit(l,n), n = 8,14)
      end do 
      close(220)

! 3. initialization of tracers

      do i = 1,iip1
         do j = 1,jjp1
            do l = 1,llm

               pgcm = aps(l) + bps(l)*ps(i,j)  ! gcm pressure
               pgcm = log(pgcm)                ! for the vertical interpolation
               mmean(i,j,l) = 0.

! 3.1 vertical interpolation

               do n = 1,nspe
                  call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt)
                  vmrint(n) = vmr
                  iq = niq(n)
                  mmean(i,j,l) = mmean(i,j,l) + vmrint(n)*mmol(iq)
               end do

! 3.2 attribute mixing ratio: - all layers or only thermosphere
!                             - with our without h2o 

               if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 1.e-3)) then
                  do n = 1,nspe
                     iq = niq(n)
                     if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
                        pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
                     end if
                  end do
               end if

            end do
         end do
      end do

      ! set surface values of chemistry tracers to zero
      if (flagthermo == 0) then
        ! NB: no problem for "surface water vapour" tracer which is always 0
        do n=1,nspe
          iq=niq(n)
          qsurf(1:ngridmx,iq)=0
        enddo
      endif


! 3.3 initialization of tracers not contained in the initialization files

! methane : 10 ppbv

      if (igcm_ch4 /= 0) then
         vmr = 10.e-9       
         do i = 1,iip1
            do j = 1,jjp1
               do l = 1,llm
                  pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l)
               end do
            end do
         end do
         ! set surface value to zero
         qsurf(1:ngridmx,igcm_ch4)=0
      end if

      end