source: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/inichim_1D.F90 @ 3233

      subroutine inichim_1D(nlayer, nq, pq, qsurf, play, &
                                  flagh2o,flagthermo)

      use tracer_h, only: noms, mmol
      use datafile_mod, only: datadir

      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
!  Rewritten 03/2021 Yassin Jaziri (Use of #Moderntrac-v1 to init thanks traceur.def)
!
!  Arguments:
!  ----------
!
!  nlayer          Number of atmospheric layers
!  pq(nlayer,nq)   Advected fields, ie chemical species here
!  qsurf(nq)       Amount of tracer on the surface (kg/m2)
!  play(nlayer)    Pressure (Pa)
!  flagh2o         flag for initialisation of h2o (1: yes / 0: no)
!  flagthermo      flag for initialisation of thermosphere only (1: yes / 0: no)
!
!=======================================================================


! inputs :

      integer,intent(in) :: nlayer         ! Number of atmospheric layers.
      integer,intent(in) :: nq             ! number of tracers
      real   ,intent(in) :: play(nlayer)   ! Mid-layer pressure (Pa).  
      integer,intent(in) :: flagh2o        ! flag for h2o initialisation
      integer,intent(in) :: flagthermo     ! flag for thermosphere initialisation only

! outputs :

      real,intent(out)   :: pq(nlayer,nq)  ! tracers (kg/kg_of_air)
      real,intent(out)   :: qsurf(nq)      ! surface values (kg/m2) of tracers

! local :

      real, allocatable  :: pf(:)          ! pressure in vmr profile files set in traceur.def
      real, allocatable  :: qf(:)          ! vmr      in vmr profile files set in traceur.def

      real    :: mmean(nlayer)             ! mean molecular mass (g)
      real    :: pqx(nlayer,nq)            ! tracers (vmr)
      real    :: qx(nq)                    ! constant vmr set in traceur.def
      integer :: iq, ilay, iline, nlines, ierr

      CHARACTER(len=100) :: qxf(nq)        ! vmr profile files set in traceur.def
      CHARACTER(len=100) :: fil            ! path files
      character(len=500) :: tracline       ! to store a line of text

      logical :: foundback = .false.

! 1. initialization

      pq(:,:)  = 0.
      qsurf(:) = 0.
      pqx(:,:) = 0.
      qx(:)    = 0.
      qxf(:)   = 'None'
      nlines   = 0

! 2. load in traceur.def chemistry data for initialization:

      ! Skip nq
      open(90,file='traceur.def',status='old',form='formatted',iostat=ierr)
      if (ierr.eq.0) then
        READ(90,'(A)') tracline
        IF (trim(tracline).eq.'#ModernTrac-v1') THEN ! Test modern traceur.def
           DO
             READ(90,'(A)',iostat=ierr) tracline
             IF (ierr.eq.0) THEN
               IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
                 EXIT
               ENDIF
             ELSE ! If pb, or if reached EOF without having found number of tracer
                write(*,*) "calchim: error reading line of tracers"
                write(*,*) "   (first lines of traceur.def) "
                stop
             ENDIF
           ENDDO
        ENDIF ! if modern or standard traceur.def
      else
         write(*,*) "calchim: error opening traceur.def in inichim_1D"
         stop
      endif

      ! Get data of tracers
      do iq=1,nq
         read(90,'(A)') tracline
         write(*,*)"inichim_1D: iq=",iq,"noms(iq)=",trim(noms(iq))
         if (index(tracline,'qx='    ) /= 0) then
            read(tracline(index(tracline,'qx=')+len('qx='):),*) qx(iq)
            write(*,*) ' Parameter value (traceur.def) : qx=', qx(iq)
         else
            write(*,*) ' Parameter value (default)     : qx=', qx(iq)
         end if
         if (index(tracline,'qxf='    ) /= 0) then
            read(tracline(index(tracline,'qxf=')+len('qxf='):),*) qxf(iq)
            write(*,*) ' Parameter value (traceur.def) : qxf=', qxf(iq)
         else
            write(*,*) ' Parameter value (default)     : qxf=', qxf(iq)
         end if
      end do

      close(90)

! 3. initialization of tracers

! 3.1 vertical interpolation

      do iq=1,nq
         if (qx(iq) /= 0.) then
            pqx(:,iq) = qx(iq)
         else if (qxf(iq) /= 'None') then
            ! Opening file
            fil = trim(datadir)//'/chemical_profiles/'//qxf(iq)
            print*, 'chemical pofile '//trim(noms(iq))//': ', fil
            open(UNIT=90,FILE=fil,STATUS='old',iostat=ierr)
            if (ierr.eq.0) then
               read(90,*) ! read one header line
               do         ! get number of lines
                   read(90,*,iostat=ierr)
                   if (ierr<0) exit
                   nlines = nlines + 1
               end do
               ! allocate reading variable
               allocate(pf(nlines))
               allocate(qf(nlines))
               ! read file
               rewind(90) ! restart from the beggining of the file
               read(90,*) ! read one header line
               do iline=1,nlines
                  read(90,*) pf(iline), qf(iline) ! pf [Pa], qf [vmr]
               end do
               ! interp in gcm grid
               do ilay=1,nlayer
                  call intrplf(log(play(ilay)),pqx(ilay,iq),log(pf),qf,nlines)
               end do
               ! deallocate for next tracer
               deallocate(pf)
               deallocate(qf)
            else
               write(*,*) 'inichim_1D: error opening ', fil
               stop
            endif
            close(90)
         end if
      end do

! 3.2 background gas

      do iq=1,nq
         if (all(pqx(:,iq)==1.)) then
            pqx(:,iq) = 0.
            do ilay=1,nlayer
               pqx(ilay,iq) = 1-sum(pqx(ilay,:))
               if (pqx(ilay,iq)<=0.) then
                  write(*,*) 'inichim_1D: vmr tot > 1 not possible'
                  stop
               end if
            end do
            foundback = .true.
            exit ! you found the background gas you can skip others
         end if
      end do
      if (.not.foundback) then
         write(*,*) 'inichim_1D: you need to set a background gas'
         write(*,*) '            by qx=1. in traceur.def'
         stop
      end if

! 3.3 mmean
      mmean(:) = 0.
      do ilay=1,nlayer
         do iq=1,nq
            mmean(ilay) = mmean(ilay) + pqx(ilay,iq)*mmol(iq)
         end do
      end do

! 3.4 convert vmr to mmr

      do ilay=1,nlayer
         do iq=1,nq
            pq(ilay,iq) = pqx(ilay,iq)*mmol(iq)/mmean(ilay)
         end do
      end do

! 4. Hard coding
! Do whatever you want here to specify pq and qsurf
! Or use #ModernTrac-v1 and add another option section 2.

      end