ecearth3/sources/tm5mp/proj/cbm4/m7_coaset.F90

#include "tm5.inc"

SUBROUTINE m7_coaset(kproma, kbdim, klev,  paernl, ptp1, &
                     papp1,  pm6rp, prhop, pcom          )
  !
  ! *m7_coaset*  calculates the coagulation kernels between the modes
  !
  ! Authors:
  ! ---------
  ! J. Wilson  and E. Vignati, JRC/EI (original source)                09/2000
  ! P. Stier, MPI                     (f90-version, changes, comments)    2001 
  !
  ! Modifications:
  ! --------------
  ! Philip Stier, MPI                             2001
  !
  ! Purpose
  ! ---------
  ! This routine calculates the coaglation kernels between particles
  ! with the count median radii of the three modes.
  ! Coagulation allowed between the following modes:
  ! soluble modes:   1+1=1, 2+2=2, 1+2=2, 1+3=3, 1+4=4, 2+3=3, 2+4=4
  ! insoluble modes: 2i+2i=2i
  ! mixed modes:     1+2i=2, 1+4i=4, 2+2i=2, 2+4i=4, 3+2i=3, 4+2i=4
  !
  ! Interface:
  ! -----------
  ! *m7_coaset* is called from *m7_dnum*
  !
  ! Externals:
  ! -----------
  ! none
  !
  ! Reference:
  ! -----------
  ! The calculations are based on:
  ! Fuchs, N.A. (1964). The Mechanics of Aerosols. Pergamon Press. Oxford. 
  ! (Chapter VII, Section 49)
  !
  !  Warning:
  !  --------
  !  For optimization purposes currently only "physically reasonable" elements of the
  !  coagulation kernel pcom are calculated in m7_concoag. These elements are specified
  !  in the matrix locoagmask in mo_aero_m7. Check carefully and adapt locoagmask 
  !  accordingly  before changes in the code below.
  !
  USE mo_aero_m7, ONLY: pi, bk, sqrt2, nmod, locoagmask
  !
  IMPLICIT NONE
  !
  !--- Parameter list:
  !
  !  paernl            = aerosol number for each mode [cm-3]
  !  ptp1              = atmospheric temperature at time t+1 [K]
  !  papp1             = atmospheric pressure at time t+1 [Pa]
  !  pm6rp             = mean mode actual radius (wet radius for soluble modes 
  !                      and dry radius for insoluble modes) [cm]
  !  prhop             = mean mode particle density [g cm-3]
  !  pcom(:,:,jm1,jm2) = Coagulation coefficient for modes jm1 and jm2 []
  !
  !--- List of local variables:
  !
  ! zwlc              = mean free pathlength []
  ! zairvisc          = air viscosity []
  ! zrpav             = average radius of the two interacting modes
  ! zpvx              = volume of the xth interacting mode
  ! zpmx              = mass of the xth interacting mode
  !
  ! zrknudx           = knudsen number of the xth interacting mode
  ! zpd2x             = particle diffusion of the xth interacting mode
  !
  !--- Parameters:
  !
  INTEGER :: kproma, kbdim, klev

  REAL    :: ptp1(kbdim,klev),          papp1(kbdim,klev)

  REAL    :: pm6rp(kbdim,klev,nmod),    prhop(kbdim,klev,nmod),                 &
             paernl(kbdim,klev,nmod)

  REAL    :: pcom(kbdim,klev,nmod,nmod)

  !--- Local variables:
  
  INTEGER :: jm2, jm1, jl, jk
  !
  REAL    :: zpbyone,     zwlc,        zairvisc,    zbtketc,   &
             zrpvm1,      zrpvm2,      zrpav,       zpv1,      &
             zpm1,        zcv21,       zpv2,        zpm2,      &
             zcv22,       zcv2av,      zrknud1,     zrknud2,   &
             ze1,         ze2,         zpd21,	    zpd22,     &
             zpdav,       zxd1,        zh21,        zxd2,      &
             zh22,        zh2av,       zcoc,        zhu1,      &
             zhu2,        zhu

  !--- 1) Calculation of the coagulation coefficient: ---------------------------
  !
  DO jm2=1,nmod
     DO jm1=jm2,nmod

        IF (locoagmask(jm1,jm2)) THEN

        DO jk=1,klev
           DO jl=1,kproma
	      IF (paernl(jl,jk,jm1) > 1.E-10 .AND. paernl(jl,jk,jm2) > 1.E-10 .AND. &
                  pm6rp(jl,jk,jm1)  > 1.E-10 .AND. pm6rp(jl,jk,jm2)  > 1.E-10 .AND. &
                  prhop(jl,jk,jm1)  > 1.E-10 .AND. prhop(jl,jk,jm2)  > 1.E-10        ) THEN        
           
                 !--- 1.1) Calculate ambient properties:

                 !--- Mean free pathlength ? (from Knudsen Number below):
                 !    Parametrisation?
                 zpbyone=1000.0 / (papp1(jl,jk)/100.0)
                 zwlc=6.6e-6 * ptp1(jl,jk) / 293.15 * zpbyone

                 !--- Viscosity:
                 zairvisc=1.827e-4 * (ptp1(jl,jk) / 291.15)**0.74

                 !--- 
                 zbtketc=bk * ptp1(jl,jk) / 6.0 / pi / zairvisc

                 !--- Count median radii of the respective modes:
                 zrpvm1=pm6rp(jl,jk,jm1)
                 zrpvm2=pm6rp(jl,jk,jm2)

                 !--- Average radius of the modes:
                 zrpav=(zrpvm1 + zrpvm2) / 2.0

                 !--- Volume and mass of mode jm1:
                 zpv1=4.0 / 3.0 * pi * zrpvm1**3.0
                 zpm1=zpv1 * prhop(jl,jk,jm1)

                 !--- Squared mean particle velocity of mode jm1:
                 zcv21=8.0 * bk * ptp1(jl,jk) / (pi * zpm1)

                 !--- Volume and mass of particles in mode jm2:
                 zpv2=4.0 / 3.0 * pi * zrpvm2**3.0
                 zpm2=zpv2 * prhop(jl,jk,jm2)

                 !--- Squared mean particle velocity of mode jm2:
                 zcv22=8.0 * bk * ptp1(jl,jk) / (pi * zpm2)

                 !--- Fuchs: G_r (below Eq. 49.27):
                 zcv2av=SQRT(zcv21 + zcv22)

                 !--- Knudsen numbers of the modes:
                 !@@@ Check: is zwlc mean free path then zrknud=zwlc/zrpvm!
                 zrknud1=zwlc/zrpvm1/2.0
                 zrknud2=zwlc/zrpvm2/2.0

                 !--- Diffusivities of the respective modes:
                 !    !@@@ Parameterisation?
                 ze1=EXP(-0.43/zrknud1)
                 ze2=EXP(-0.43/zrknud2)
                 zpd21=zbtketc * (1.0 + zrknud1*2.492 + zrknud1*0.84*ze1) / zrpvm1
                 zpd22=zbtketc * (1.0 + zrknud2*2.492 + zrknud2*0.84*ze2) / zrpvm2

                 !--- Average diffusivity of the modes:
                 zpdav=(zpd21 + zpd22) / 2.0

                 !--- Average mean free path of particles in jm1:
                 zxd1=8.0 * zpd21 / (pi*SQRT(zcv21))

                 !--- Mean distance from surface after mean free path (Eq. 49.13):
                 zh21=(((2.0*zrpvm1 + zxd1)**3.0 -                        &
                       (4.0*zrpvm1*zrpvm1 + zxd1*zxd1)**1.5) /            &
                       (6.0*zrpvm1*zxd1) - 2*zrpvm1           ) * sqrt2

                 !--- Average mean free path of particles in jm2:
                 zxd2=8.0 * zpd22 / (pi*SQRT(zcv22))

                 !--- Mean distance from surface after mean free path (Eq. 49.13):

                 zh22=(((2.0*zrpvm2 + zxd2)**3.0 -                        &
                       (4.0*zrpvm2*zrpvm2 + zxd2*zxd2)**1.5) /            &
                       (6.0*zrpvm2*zxd2) - 2*zrpvm2           ) * sqrt2

                 !--- Fuchs: delta_r !@@@ (why division by sqrt2?)
                 zh2av=SQRT(zh21*zh21 + zh22*zh22) / sqrt2
                 
                 !--- 1.2) Calculation of the coagulation coefficient pcom (Eq. 49.26):
                 !         Factor 16 instead factor 8 as in Fuchs as his formulation
                 !         applies for the inter-modal coagulation. This is taken into
                 !         account in the assignment of the inter-modal coagulation
                 !         coefficient.

                 zcoc=16.0 * pi * zpdav * zrpav

                 !--- Calculation of beta=1/zhu (Eq. 49.27):
                 zhu1=4.0 * zpdav / (zcv2av * zrpav)
                 zhu2=zrpav / (zrpav + zh2av /2.0)
                 zhu=zhu1 +  zhu2

                 !--- Coagulation coefficient following (Eq.49.26):
                 pcom(jl,jk,jm1,jm2)=zcoc / zhu

	      ELSE
                 pcom(jl,jk,jm1,jm2)=0.
	      END IF

              !--- 2) Mirror the coagulation matrix (symmetric): -----------------

              pcom(jl,jk,jm2,jm1)=pcom(jl,jk,jm1,jm2)

           END DO
        END DO

        ELSE

           pcom(1:kproma,:,jm1,jm2)=0.

        END IF ! locoagmask

     END DO
  END DO

END SUBROUTINE m7_coaset