source: trunk/LMDZ.PLUTO/libf/muphypluto/mp2m_methods.F90 @ 3607

MODULE MP2M_METHODS
    !============================================================================
    !
    !     Purpose
    !     -------
    !     Model miscellaneous methods module.
    !
    !     The module contains miscellaneous methods used in the haze of the model.
    !     The module contains nine methods:
    !        - mm_lambda_air
    !        - mm_eta_air
    !        - mm_ps2s
    !        - mm_qmean
    !        - mm_get_btk
    !        - mm_get_kco
    !        - mm_get_kfm
    !
    !     Authors
    !     -------
    !     B. de Batz de Trenquelléon, J. Burgalat (11/2024)
    !
    !============================================================================

    USE MP2M_MPREC
    USE MP2M_GLOBALS
    USE LINT_DSET
    USE LINT_LOCATORS
    IMPLICIT NONE

    PRIVATE

    PUBLIC  :: mm_lambda_air, mm_eta_air, mm_ps2s, mm_qmean, mm_get_btk, mm_get_kfm, mm_get_kco


    CONTAINS

    !============================================================================
    !                          GENERAL METHODS
    !============================================================================

    ELEMENTAL FUNCTION mm_lambda_air(T,P) RESULT(res)
      !! Get the air mean free path at given temperature and pressure.
      !!
      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
      REAL(kind=mm_wp), INTENT(in) :: P ! Pressure level (Pa).
      REAL(kind=mm_wp) :: res           ! Air mean free path (m).

      res = (mm_kboltz*T) / (dsqrt(2._mm_wp)*mm_pi*(2._mm_wp*mm_air_rad)**2*P)
      
      RETURN
    END FUNCTION mm_lambda_air

    
    ELEMENTAL FUNCTION mm_eta_air(T) RESULT (res)
      !! Get the air dynamical viscosity at a given temperature using Sutherland method.
      !!
      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
      REAL(kind=mm_wp) :: res           ! Air viscosity at given temperature (Pa.s-1).

      REAL(kind=mm_wp), PARAMETER :: eta0 = 1.74e-5_mm_wp
      REAL(kind=mm_wp), PARAMETER :: Tsut = 109._mm_wp
      REAL(kind=mm_wp), PARAMETER :: Tref = 293._mm_wp
      
      res = eta0 * dsqrt(T/Tref) * ((1._mm_wp + Tsut/Tref) / (1._mm_wp + Tsut/T))
      
      RETURN
    END FUNCTION mm_eta_air


    !============================================================================
    !                      AEROSOL COAGULATION METHODS
    !============================================================================
    
    FUNCTION mm_ps2s(rcs,k,flow) RESULT(res)
      !! Get the proportion of aerosols that remains in the spherical mode during SS coagulation.
      !!
      !! From __k__ and __flow__ values, the method selects one of the four probability datasets
      !! in mm_globals(module) module (for instance mm_pco0p) and interpolates linearly probability
      !! for the given value of __rcs__, __T__ and __P__.
      !!
      !! @Warning
      !! Here, the method assumes the datasets define the probability for __spherical__ particles to
      !! be transferred in the __fractal__ mode, but returns the proportion of particles that remains
      !! in the mode (which is expected by MP2M model).
      !!
      ! Characteristic radius of the spherical size-distribution (m).
      REAL(kind=mm_wp), INTENT(in) :: rcs
      ! Order of the moment (0 or 3 expected).
      INTEGER, INTENT(in)          :: k
      ! Flow regime indicator (0: Continuous - Kn << 1, 1: Free-Molecular - Kn >> 1).
      INTEGER, INTENT(in)          :: flow
      
      ! Proportion of spherical particles that stay in the spherical mode during SS coagulation.
      REAL(kind=mm_wp) :: res
      
      ! Local variable.
      TYPE(dset1d), POINTER :: pp
      
      res = 1._mm_wp
      IF (rcs <= 0.0_mm_wp .OR. .NOT.mm_w_ps2s) RETURN
      
      SELECT CASE(k+flow)
        CASE(0)      ; pp => mm_pco0p ! 0 = 0 + 0 -> M0 / CO
        CASE(1)      ; pp => mm_pfm0p ! 1 = 0 + 1 -> M0 / FM
        CASE(3)      ; pp => mm_pco3p ! 3 = 3 + 0 -> M3 / CO
        CASE(4)      ; pp => mm_pfm3p ! 4 = 3 + 1 -> M3 / FM
        CASE DEFAULT ; RETURN
      END SELECT

      IF (.NOT.hdcd_lint_dset(rcs,pp,locate_reg_ext,res)) THEN
        WRITE(*,'(a)') "mp2m_methods:ps2s_sc: Cannot interpolate transfert probability"
        call EXIT(10)
      ELSE
        ! Sanity check: bound probability value between 0 and 1.
        res = MAX(0.0_mm_wp,MIN(res,1.0_mm_wp))
        ! We have interpolated f = 1 - p and we need p !
        res = 1._mm_wp - res
      ENDIF
    END FUNCTION mm_ps2s


    FUNCTION mm_qmean(rc1,rc2,order,modes,T) RESULT(res)
      !! Get the electric correction for coagulation kernel.
      !!
      !! The method computes the eletric charging correction to apply to the coagulation
      !! kernel as a function of the temperature and the characteristic radius of the
      !! mode involved in the coagulation.
      !! Here the electric charging correction is computed using linear interpolation from
      !! pre-tabulated values.
      !!
      !! @Warning:
      !! Modes are referred by a two letters uppercase string with the combination of:
      !!    - S : spherical mode
      !!    - F : fractal mode
      !!
      REAL(kind=mm_wp), INTENT(in) :: rc1   ! Characteristic radius of the first mode (m).
      REAL(kind=mm_wp), INTENT(in) :: rc2   ! Characteristic radius of the the second mode (m).
      INTEGER, INTENT(in)          :: order ! Moment's order (0 or 3 expected).
      CHARACTER(len=2), INTENT(in) :: modes ! Interaction mode (combination of [S,F]).
      REAL(kind=mm_wp), INTENT(in) :: T     ! Temperature (K).
      
      ! Electric charging correction.
      REAL(kind=mm_wp) :: res

      ! Local variable.
      INTEGER          :: chx
      REAL(kind=mm_wp) :: r_tmp, t_tmp
      
      chx = 0
      IF (.NOT.mm_w_qe) THEN
        res = 1._mm_wp
        RETURN
      ENDIF
    
      IF (SCAN(modes(1:1),"sS") /= 0) chx = chx + 1
      IF (SCAN(modes(2:2),"sS") /= 0) chx = chx + 1
      IF (SCAN(modes(1:1),"fF") /= 0) chx = chx + 3
      IF (SCAN(modes(2:2),"fF") /= 0) chx = chx + 3
      
      chx = chx + order
      
      SELECT CASE(chx)
        ! M0/SS:
        CASE(2)
          res = 1._mm_wp
        ! M0/SF:
        CASE(4)     
          ! Fix max values of input parameters
          r_tmp = MAX(MIN(log(rc1),mm_qbsf0_e(2,2)),mm_qbsf0_e(2,1))
          t_tmp = MAX(MIN(T,mm_qbsf0_e(1,2)),mm_qbsf0_e(1,1))
          ! Interpolates values
          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbsf0,locate_reg,res)) THEN
            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
            call EXIT(10)
          ENDIF
        ! M3/SS:
        CASE(5)
          res = 1._mm_wp
        ! M0/FF:
        CASE(6)
          ! Fix max values of input parameters
          r_tmp = MAX(MIN(log(rc1),mm_qbff0_e(2,2)),mm_qbff0_e(2,1))
          t_tmp = MAX(MIN(T,mm_qbff0_e(1,2)),mm_qbff0_e(1,1))
          ! Interpolates values
          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbff0,locate_reg,res)) THEN
            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
            call EXIT(10)
          ENDIF
        ! M3/SF:
        CASE(7)
          ! Fix max values of input parameters
          r_tmp = MAX(MIN(log(rc1),mm_qbsf3_e(2,2)),mm_qbsf3_e(2,1))
          t_tmp = MAX(MIN(T,mm_qbsf3_e(1,2)),mm_qbsf3_e(1,1))
          ! Interpolates values
          IF (.NOT.hdcd_lint_dset(t_tmp,r_tmp,mm_qbsf3,locate_reg,res)) THEN
            WRITE(*,'(a)') "mp2m_methods:mm_qmean: Cannot interpolate mean Qelec"
            call EXIT(10)
          ENDIF
        ! Anything else:
        CASE DEFAULT
          res = 1._mm_wp
      END SELECT
      
      RETURN
    END FUNCTION mm_qmean


    PURE FUNCTION mm_get_btk(t,k) RESULT(res)
      !! Get the value of the Free-Molecular regime coagulation pre-factor b^t_k.
      !!
      !! @Note
      !! __k__ can only be one of the following value: 0 or 3. __t__ ranges only from 1 to 5.
      !!
      INTEGER, INTENT(in) :: t ! Index of the b^t_k coefficient.
      INTEGER, INTENT(in) :: k ! Moment order of the b^t_k coefficient.
      
      ! b^t_k coefficient.
      REAL(kind=mm_wp) :: res

      ! Sanity check:
      IF (.NOT.(k == 3 .OR. k == 0)) THEN
        res = 0._mm_wp
      ENDIF
      
      ! Sanity check:
      IF (t > 5 .OR. t < 1) THEN
        res = 0._mm_wp
      ENDIF
      
      IF (k == 0) THEN
        res = mm_bt0(t)

      ELSE IF (k == 3) THEN
        res = mm_bt3(t)
      ENDIF
      
      RETURN
    END FUNCTION mm_get_btk


    ELEMENTAL FUNCTION mm_get_kco(T) RESULT(res)
      !! Get the Continuous regime (Kn << 1) thermodynamics pre-factor of the coagulation kernel.
      !!
      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
      REAL(kind=mm_wp) :: res           ! Continuous regime thermodynamics pre-factor (m3.s-1).

      res = (2._mm_wp*mm_kboltz*T) / (3._mm_wp*mm_eta_air(T))
      
      RETURN
    END FUNCTION mm_get_kco


    ELEMENTAL FUNCTION mm_get_kfm(T) RESULT(res)
      !! Get the Free-Molecular regime (Kn >> 1) thermodynamics pre-factor of the coagulation kernel.
      !!
      REAL(kind=mm_wp), INTENT(in) :: T ! Temperature (K).
      REAL(kind=mm_wp) :: res           ! Free-Molecular regime thermodynamics pre-factor (m^(5/2).s-1).

      res = (6._mm_wp*mm_kboltz*T / mm_rhoaer)**(0.5_mm_wp)
      
      RETURN
    END FUNCTION mm_get_kfm

END MODULE MP2M_METHODS