ecearth3/sources/tm5mp/proj/cb05/m7_dgas_org.F90

#include "tm5.inc"

#ifdef with_budgets
SUBROUTINE m7_dgas_org(kproma, kbdim, klev,  pelvoc, psvoc,  paerml, paernl, &
                   ptp1,   papp1, pm6rp, &
                   ptime,pprocess                                 )
#else
SUBROUTINE m7_dgas_org(kproma, kbdim, klev,  pelvoc, psvoc,  paerml, paernl, &
                   ptp1,   papp1, pm6rp, &
                   ptime                                           )
#endif
  !
  !**** *m7_dgas_org*  calculates the transfer of mass due to 
  !                organic condensation
  !
  !    Authors:
  !    -----------
  !    J. Wilson, E. Vignati, JRC/EI (original source)                05/2000
  !    R. Makkonen, UHEL                                              2015
  !

  USE mo_time_control,  ONLY: delta_time
  USE mo_control,       ONLY: nrow
  USE mo_aero_m7,       ONLY: pi,      wh2so4,  rerg,    avo,            &
                              fmax,    caccso4,                          &
                              nmod,    naermod,                          &
                              iocks,  iocas,  ioccs,  iocki,     &
                              isoans, isoaks, isoaas, isoacs, isoaki
  USE mo_aero_m7,       ONLY: condensation_sink

!  USE mo_aero_mem,      ONLY: d_cond_so4

#ifdef with_budgets
  Use M7_Data, only: nm7procs
#endif

  IMPLICIT NONE
  !
  !--- Parameter list:
  !
  ! pso4g     = mass of gas phase sulfate [molec. cm-3]
  ! pm6rp     = mean mode actual radius (wet radius for soluble modes 
  !             and dry radius for insoluble modes) [cm]
  ! pso4_x    = mass of sulphate condensed on insoluble mode x [molec. cm-3]
  ! 
  !--- Local Variables:
  !
  ! zde       = molecular diffusion []
  ! zvelb     = velocity []
  ! zcondo    = condensation coefficient []
  ! zc2(nmod) = flux of sulfate condensing on the respective mode 
  !             per sulfate gas phase concentration []
  ! zcondo    = total flux of condensing sulfate 
  !             per sulfate gas phase concentration []
  ! zfcond    = total mass of condensing sulfate for one timestep []

  INTEGER :: kproma, kbdim, klev

  REAL    :: ptime, time_step_len

  REAL    :: ptp1(kbdim,klev),        papp1(kbdim,klev),                 &
             pelvoc(kbdim,klev),      psvoc(kbdim,klev) 
  
  REAL    :: paernl(kbdim,klev,nmod), paerml(kbdim,klev,naermod),        &
             pm6rp(kbdim,klev,nmod)

#ifdef with_budgets
  Real :: pprocess(kbdim,klev,nm7procs)
#endif
  !
  ! Local variables:

  INTEGER :: jl,         jk,         jmod,       jrow
  
  REAL    :: zfcond,     zftot,      zpbyone,    zde2,                   & 
             zvelb,      zxibc,      zm6rp,      zf1,                    &
             zqtmst

  REAL    :: zcondo(kbdim,klev)

  REAL    :: zc2(kbdim,klev,nmod)

  REAL    :: lambda, cs(7), cs_sum, modeweight(7), modeweight_sum
!  REAL(dp):: soa_yield, org, lambda, cs(7), cs_weighted(7), cssum, pblheight, oc_mass_sum
  REAL, parameter  ::  oc2pom_soa = 2.4   !RM test 1e1

  !--- 0) Initialisations: -------------------------------------------------
  !
  jrow=nrow(2)

  zcondo(:,:)=0.0
      
  zc2(:,:,:) = 0.0  

  time_step_len = ptime

  zqtmst=1/time_step_len
  

  !--- 1) Calculate condensation rate for cm diameter sulphate aerosols: ---
  !
  DO jmod=1,nmod
     DO jk=1,klev
        DO jl=1,kproma
           IF (pm6rp(jl,jk,jmod).GT.0.) THEN

              !--- Diffusion coefficient (Reference???):

              zpbyone=1000.0 / (papp1(jl,jk)/100.0)

              zde2=0.073 * zpbyone * (ptp1(jl,jk) / 298.15)**1.5  

              !--- Mean molecule velocity (Moore, 1962 (S+P equ. 8.2)):

              zvelb=SQRT(8.0 * rerg * ptp1(jl,jk) / pi / wh2so4)  

              !--- ???Fuchs???

              zxibc=8.0 * zde2 / pi / zvelb
              !
              ! Use count median radius:

              zm6rp=pm6rp(jl,jk,jmod)

              !--- Distance from particle up to which the kinetic regime applies:

              zf1=( (zm6rp + zxibc)**3.0 - (zm6rp**2.0 + zxibc**2.0)**1.5 ) / &
                  (3.0 * zm6rp * zxibc) - zm6rp

              !--- Diffusive flux to single particle surface:
              !    (Elisabetta's thesis: fraction in equ. 2.26)

              zc2(jl,jk,jmod)=(4.0 * pi * zde2 * zm6rp ) /                      &
                              ((4.0 * zde2) / (zvelb * zm6rp * caccso4(jmod)) +   &
                               (zm6rp/(zm6rp+zf1))                              )

              !--- Total diffusive flux to all particles in the respective mode:
              !    (per concentration of gas phase sulfate)

              zc2(jl,jk,jmod)=zc2(jl,jk,jmod) * paernl(jl,jk,jmod)

              !--- Total diffusive flux to all particles in all modes:
              !    (per concentration of gas phase sulfate)

              zcondo(jl,jk)=zcondo(jl,jk)+ zc2(jl,jk,jmod)  

           END IF
        END DO
     END DO
  END DO
  !
  !--- 2) Calculation of the new organic aerosol masses and of the ---------
  !       mass of organic condensing on the respective modes:
  !
  DO jk=1,klev
     DO jl=1,kproma

        lambda=6.6E-8

        !
        !--- Calculate condensation sink and OA-mass of each mode
        !

        CALL condensation_sink(7,paernl(jl,jk,:), 0.01*pm6rp(jl,jk,:), lambda, cs)

        modeweight(1)=0.
        modeweight(2)=MIN(MAX(paerml(jl,jk,iocks)+paerml(jl,jk,isoaks),0.),1.E5)
        modeweight(3)=MIN(MAX(paerml(jl,jk,iocas)+paerml(jl,jk,isoaas),0.),1.E5)
        modeweight(4)=MIN(MAX(paerml(jl,jk,ioccs)+paerml(jl,jk,isoacs),0.),1.E5)
        modeweight(5)=MIN(MAX(paerml(jl,jk,iocki)+paerml(jl,jk,isoaki),0.),1.E5)
        modeweight(6)=0.
        modeweight(7)=0.

        cs_sum=SUM(cs(1:5))  ! Sum only over OC-modes
        modeweight_sum=SUM(modeweight)

        !
        !--- Condense S/LVOC according to mode OA mass
        !

        IF(modeweight_sum .GT. 1.E-13 .AND. psvoc(jl,jk).GT.1.e-16 .AND. psvoc(jl,jk).LT.1.e10) THEN

          paerml(jl,jk,isoaks)  = paerml(jl,jk,isoaks)  + psvoc(jl,jk)*modeweight(2)/modeweight_sum
          paerml(jl,jk,isoaas)  = paerml(jl,jk,isoaas)  + psvoc(jl,jk)*modeweight(3)/modeweight_sum
          paerml(jl,jk,isoacs)  = paerml(jl,jk,isoacs)  + psvoc(jl,jk)*modeweight(4)/modeweight_sum
          paerml(jl,jk,isoaki)  = paerml(jl,jk,isoaki)  + psvoc(jl,jk)*modeweight(5)/modeweight_sum

          pprocess(jl,jk,80)  = psvoc(jl,jk)*modeweight(2)/modeweight_sum
          pprocess(jl,jk,81)  = psvoc(jl,jk)*modeweight(3)/modeweight_sum
          pprocess(jl,jk,82)  = psvoc(jl,jk)*modeweight(4)/modeweight_sum
          pprocess(jl,jk,83)  = psvoc(jl,jk)*modeweight(5)/modeweight_sum

        !--- If only small amount of existing OA, condensing to accumulation mode

        ELSE IF(psvoc(jl,jk).GT.1.e-16 .AND. psvoc(jl,jk).LT.1.e10) THEN

          paerml(jl,jk,isoaas)  = paerml(jl,jk,isoaas)  + psvoc(jl,jk)

          pprocess(jl,jk,81)    = psvoc(jl,jk) ! Condensation 3 SOA

        END IF

        !
        !--- Condense ELVOC according to mode condensation sink
        !

        IF(cs_sum.GT.1.E-15 .AND. cs_sum.LT.1.E5 .AND. pelvoc(jl,jk).GT.1.e-16 .AND. pelvoc(jl,jk).LT.1.e10) THEN

           paerml(jl,jk,isoans)=paerml(jl,jk,isoans)+ &
                                pelvoc(jl,jk)*cs(1)/cs_sum
           paerml(jl,jk,isoaks)=paerml(jl,jk,isoaks)+ &
                                pelvoc(jl,jk)*cs(2)/cs_sum
           paerml(jl,jk,isoaas)=paerml(jl,jk,isoaas)+ &
                                pelvoc(jl,jk)*cs(3)/cs_sum
           paerml(jl,jk,isoacs)=paerml(jl,jk,isoacs)+ &
                                pelvoc(jl,jk)*cs(4)/cs_sum
           paerml(jl,jk,isoaki)=paerml(jl,jk,isoaki)+ &
                                pelvoc(jl,jk)*cs(5)/cs_sum

           pprocess(jl,jk,79) = pelvoc(jl,jk)*cs(1)/cs_sum ! Condensation 1 SOA
           pprocess(jl,jk,80) = pelvoc(jl,jk)*cs(2)/cs_sum ! Condensation 2 SOA
           pprocess(jl,jk,81) = pelvoc(jl,jk)*cs(3)/cs_sum ! Condensation 3 SOA
           pprocess(jl,jk,82) = pelvoc(jl,jk)*cs(4)/cs_sum ! Condensation 4 SOA
           pprocess(jl,jk,83) = pelvoc(jl,jk)*cs(5)/cs_sum ! Condensation 5 SOA
!           write(*,*) 'NUC AIT ACC COA AITI',zc2(jl,jk,iocks)/zcondo(jl,jk),zc2(jl,jk,iocas)/zcondo(jl,jk),zc2(jl,jk,ioccs)/zcondo(jl,jk),zc2(jl,jk,iocki)/zcondo(jl,jk)
!           write(*,*) 'NUC AIT ACC COA AITI',cs(1)/cs_sum,cs(2)/cs_sum,cs(3)/cs_sum,cs(4)/cs_sum,cs(5)/cs_sum


        !--- If only small CS, condensing to accumulation mode

        ELSE IF(pelvoc(jl,jk).GT.1.e-16 .AND. pelvoc(jl,jk).LT.1.e10) THEN

          paerml(jl,jk,isoaas)  = paerml(jl,jk,isoaas)  + pelvoc(jl,jk)

          pprocess(jl,jk,81)    = pelvoc(jl,jk) ! Condensation 3 SOA
        END IF
     END DO
  END DO

!write(*,*) 'elvoc:',pelvoc(2100,1)
!write(*,*) 'paerml iocks:',paerml(2100,1,iocks)

!  write(*,*) 'dgas', 'so4gnew', pso4g(2100,1), 'condso4', zcondo(2100,1)*pso4g(2100,1)*time_step_len, 'limit=', pso4g(2100,1)*fmax
!  write(*,*) 'cond5= ', pso4_5(2100,1), 'rad5= ', pm6rp(2100,1,5)

END SUBROUTINE m7_dgas_org