source: trunk/LMDZ.GENERIC/libf/phystd/dyn1d/initracer_1D.F90 @ 3753

!module initracer_1D_mod

!implicit none

    subroutine initracer_1D(nlayer, nq, pq, qsurf, play)

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

    implicit none

    !=======================================================================
    !
    !  Purpose:
    !  --------
    !
    !  Initialization of the tracer values reading flags from the traceur.def
    !  
    !  Authors: 
    !  -------
    !  Adapted 04/2025 from inichim_1D.F90 by Maxime Maurice
    !
    !
    !  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)
    !
    !=======================================================================


    ! 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).  

    ! 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.

    ! initialization

    ! pq(:,:)  = 0. <- q has already been initialized
    qsurf(:) = 0.
    pqx(:,:) = 0.
    qx(:)    = 0.
    qxf(:)   = 'None'
    nlines   = 0

    ! 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 initracer_1D"
        stop
    endif

    ! Get data of tracers
    do iq=1,nq
        read(90,'(A)') tracline
        write(*,*)"initracer_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)

    ! initialization of tracers

    ! vertical interpolation

    do iq=1,nq
        if (qx(iq) /= 0.) then
            do ilay=1,nlayer
                if (pq(ilay,iq) /= 0.) then
                    print*, "initracer_1D: Tracer "//trim(noms(iq))//" seems to be already initialized"
                    print*, "              and you are now trying set it a homogeneous value (flag qx)"
                    call abort
                end if
            end do
            pqx(:,iq) = qx(iq)
        else if (qxf(iq) /= 'None') then
            do ilay=1,nlayer
                if (pq(ilay,iq) /= 0.) then
                    print*, "initracer_1D: Tracer "//trim(noms(iq))//" seems to be already initialized"
                    print*, "              and you are now trying set it a homogeneous value (flag qx)"
                    call abort
                end if
            end do
            ! Opening file
            fil = qxf(iq)
            print*, 'chemical profile '//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(*,*) 'initracer_1D: error opening ', fil
                stop
            endif
            close(90)
        end if
    end do

    ! If a background gas is defined (with qx=1)

    do iq=1,nq
        if (all(pqx(:,iq)==1.)) then
           write(*,*) 'initracer_1D: gas '//trim(noms(iq))//' defined as background'
           write(*,*) '              with qx=1. in traceur.def. Its VMR is set to'
           write(*,*) '              1-sum(VMR(others)).'
           pqx(:,iq) = 0.
           do ilay=1,nlayer
              pqx(ilay,iq) = 1-sum(pqx(ilay,:))
              if (pqx(ilay,iq)<=0.) then
                 write(*,*) 'initracer_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

    ! convert vmr to mmr

    do ilay=1,nlayer
        do iq=1,nq
            pq(ilay,iq) = pqx(ilay,iq)*mmol(iq)/mugaz
        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 subroutine initracer_1D

!end module initracer_1D_mod