Index: trunk/LMDZ.TITAN/libf/muphytitan/mm_methods.f90
===================================================================
--- trunk/LMDZ.TITAN/libf/muphytitan/mm_methods.f90	(revision 1897)
+++ trunk/LMDZ.TITAN/libf/muphytitan/mm_methods.f90	(revision 3083)
@@ -1,21 +1,21 @@
-! Copyright 2013-2015,2017 Université de Reims Champagne-Ardenne 
-! Contributor: J. Burgalat (GSMA, URCA)
+! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
+! Contributors: J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
 ! email of the author : jeremie.burgalat@univ-reims.fr
-! 
+!
 ! This software is a computer program whose purpose is to compute
 ! microphysics processes using a two-moments scheme.
-! 
+!
 ! This library is governed by the CeCILL-B license under French law and
-! abiding by the rules of distribution of free software.  You can  use, 
+! abiding by the rules of distribution of free software.  You can  use,
 ! modify and/ or redistribute the software under the terms of the CeCILL-B
 ! license as circulated by CEA, CNRS and INRIA at the following URL
-! "http://www.cecill.info". 
-! 
+! "http://www.cecill.info".
+!
 ! As a counterpart to the access to the source code and  rights to copy,
 ! modify and redistribute granted by the license, users are provided only
 ! with a limited warranty  and the software's author,  the holder of the
 ! economic rights,  and the successive licensors  have only  limited
-! liability. 
-! 
+! liability.
+!
 ! In this respect, the user's attention is drawn to the risks associated
 ! with loading,  using,  modifying and/or developing or reproducing the
@@ -25,8 +25,8 @@
 ! professionals having in-depth computer knowledge. Users are therefore
 ! encouraged to load and test the software's suitability as regards their
-! requirements in conditions enabling the security of their systems and/or 
-! data to be ensured and,  more generally, to use and operate it in the 
-! same conditions as regards security. 
-! 
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and,  more generally, to use and operate it in the
+! same conditions as regards security.
+!
 ! The fact that you are presently reading this means that you have had
 ! knowledge of the CeCILL-B license and that you accept its terms.
@@ -35,5 +35,6 @@
 !! summary: Model miscellaneous methods module.
 !! author: J. Burgalat
-!! date: 2013-2015,2017
+!! date: 2013-2015,2017,2022
+!! corrections: B. de Batz de Trenquelléon (2023)
 
 MODULE MM_METHODS
@@ -43,10 +44,10 @@
   !!
   !! All thermodynamic functions related to cloud microphysics (i.e. [[mm_methods(module):mm_lHeatX(interface)]],
-  !! [[mm_methods(module):mm_sigX(interface)]] and [[mm_methods(module):mm_psatX(interface)]]) compute related equations 
+  !! [[mm_methods(module):mm_sigX(interface)]] and [[mm_methods(module):mm_psatX(interface)]]) compute related equations
   !! from \cite{reid1986}. A version of the book is freely available [here](http://f3.tiera.ru/3/Chemistry/References/Poling%20B.E.,%20Prausnitz%20J.M.,%20O'Connell%20J.P.%20The%20Properties%20of%20Gases%20and%20Liquids%20(5ed.,%20MGH,%202000)(ISBN%200070116822)(803s).pdf).
   !!
   !! The module defines the following functions/subroutines/interfaces:
   !!
-  !! | name        | description 
+  !! | name        | description
   !! | :---------: | :-------------------------------------------------------------------------------------
   !! | mm_lheatx   | Compute latent heat released
@@ -64,8 +65,8 @@
   IMPLICIT NONE
 
-  PRIVATE 
-
-  PUBLIC  :: mm_sigX, mm_LheatX, mm_psatX, mm_qsatx, mm_fshape, &
-             mm_get_kco, mm_get_kfm, mm_eta_g, mm_lambda_g
+  PRIVATE
+
+  PUBLIC  :: mm_sigX, mm_LheatX, mm_psatX, mm_qsatx, mm_ysatX, mm_fshape, &
+    mm_get_kco, mm_get_kfm, mm_eta_g, mm_lambda_g
 
   ! ---- INTERFACES
@@ -78,13 +79,13 @@
   !! FUNCTION mm_sigX(temp,xESP)
   !! ```
-  !! 
-  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
+  !!
+  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
   !! computes the result for all the temperatures and returns a vector of same size than __temp__.
   INTERFACE mm_sigX
     MODULE PROCEDURE sigx_sc,sigx_ve
-  END INTERFACE
+  END INTERFACE mm_sigX
 
   !> Interface to Latent heat computation functions.
-  !! 
+  !!
   !! The method computes the latent heat released of a given specie at given temperature(s).
   !!
@@ -93,9 +94,9 @@
   !! ```
   !!
-  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
+  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
   !! computes the result for all the temperatures and returns a vector of same size than __temp__.
   INTERFACE mm_LheatX
     MODULE PROCEDURE lheatx_sc,lheatx_ve
-  END INTERFACE
+  END INTERFACE mm_LheatX
 
   !> Interface to saturation vapor pressure computation functions.
@@ -104,28 +105,44 @@
   !! FUNCTION mm_psatX(temp,xESP)
   !! ```
-  !! 
+  !!
   !! The method computes the saturation vapor pressure of a given specie at given temperature(s).
   !!
-  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method 
+  !! __xESP__ must always be given as a scalar. If __temp__ is given as a vector, then the method
   !! computes the result for all the temperatures and returns a vector of same size than __temp__.
   INTERFACE mm_psatX
     MODULE PROCEDURE psatx_sc,psatx_ve
-  END INTERFACE
-
-  !! Interface to saturation mass mixing ratio computaiton functions.
-  !!
-  !! The method computes the mass mixing ratio at saturation of a given specie at given temperature(s) 
+  END INTERFACE mm_psatX
+
+  !> Interface to saturation mass mixing ratio computation functions.
+  !!
+  !! The method computes the mass mixing ratio at saturation of a given specie at given temperature(s)
   !! and pressure level(s).
   !!
   !! ```fortran
-  !! FUNCTION mm_qsatX(temp,pres,xESP) 
-  !! ```
-  !!
-  !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same 
+  !! FUNCTION mm_qsatX(temp,pres,xESP)
+  !! ```
+  !!
+  !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same
   !! size !), then the method computes the result for each couple of (temperature, pressure) and returns
   !! a vector of same size than __temp__.
   INTERFACE mm_qsatx
     MODULE PROCEDURE qsatx_sc,qsatx_ve
-  END INTERFACE
+  END INTERFACE mm_qsatx
+
+  !> Interface to saturation molar mixing ratio computation functions.
+  !!
+  !! The method computes the molar mixing ratio at saturation of a given specie at given temperature(s)
+  !! and pressure level(s) [Fray and Schmidt (2009)].
+  !!
+  !! ```fortran
+  !! FUNCTION mm_ysatX(temp,pres,xESP)
+  !! ```
+  !!
+  !! __xESP__ must always be given as a scalar. If __temp__ and __pres__  are given as a vector (of same
+  !! size !), then the method computes the result for each couple of (temperature, pressure) and returns
+  !! a vector of same size than __temp__.
+  INTERFACE mm_ysatX
+    MODULE PROCEDURE ysatX_sc,ysatX_ve
+  END INTERFACE mm_ysatX
 
   !> Interface to shape factor computation functions.
@@ -137,17 +154,17 @@
   !! ```
   !!
-  !! Where __m__ is cosine of the contact angle and __x__ the curvature radius. __m__ must always be 
-  !! given as a scalar. If __x__ is given as a vector, then the method compute the result for each 
+  !! Where __m__ is cosine of the contact angle and __x__ the curvature radius. __m__ must always be
+  !! given as a scalar. If __x__ is given as a vector, then the method compute the result for each
   !! value of __x__ and and returns a vector of same size than __x__.
   INTERFACE mm_fshape
     MODULE PROCEDURE fshape_sc,fshape_ve
-  END INTERFACE
-
-  CONTAINS
+  END INTERFACE mm_fshape
+
+CONTAINS
 
   FUNCTION fshape_sc(cost,rap) RESULT(res)
     !! Get the shape factor of a ccn (scalar).
     !!
-    !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle. 
+    !! The method computes the shape factor for the heterogeneous nucleation on a fractal particle.
     !! Details about the shape factor can be found in \cite{prup1978}.
     REAL(kind=mm_wp), INTENT(in) :: cost !! Cosine of the contact angle.
@@ -160,5 +177,5 @@
       phi = dsqrt(1._mm_wp-2._mm_wp*cost*rap+rap**2)
       a = 1._mm_wp + ( (1._mm_wp-cost*rap)/phi )**3
-      b = (rap**3) * (2._mm_wp-3._mm_wp*(rap-cost)/phi+((rap-cost)/phi)**3)
+      b = (rap**3) * (2._mm_wp - 3._mm_wp*(rap-cost)/phi + ((rap-cost)/phi)**3)
       c = 3._mm_wp * cost * (rap**2) * ((rap-cost)/phi-1._mm_wp)
       res = 0.5_mm_wp*(a+b+c)
@@ -177,5 +194,5 @@
     WHERE(rap > 3000._mm_wp)
       res = ((2._mm_wp+cost)*(1._mm_wp-cost)**2)/4._mm_wp
-    ELSEWHERE 
+    ELSEWHERE
       phi = dsqrt(1._mm_wp-2._mm_wp*cost*rap+rap**2)
       a = 1._mm_wp + ((1._mm_wp-cost*rap)/phi )**3
@@ -192,13 +209,13 @@
     !! The method computes the latent heat equation as given in \cite{reid1986} p. 220 (eq. 7-9.4).
     IMPLICIT NONE
-    ! - DUMMY 
+    ! - DUMMY
     REAL(kind=mm_wp), INTENT(in) :: temp !! temperature (K).
     TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
     REAL(kind=mm_wp) :: res              !! Latent heat of given specie at given temperature (\(J.kg^{-1}\)).
     REAL(kind=mm_wp) :: ftm
-    ftm=MIN(1._mm_wp-temp/xESP%tc,1.e-3_mm_wp)
+    ftm=MAX(1._mm_wp-temp/xESP%tc,1.e-3_mm_wp)
     res = mm_rgas*xESP%tc*(7.08_mm_wp*ftm**0.354_mm_wp+10.95_mm_wp*xESP%w*ftm**0.456_mm_wp)/xESP%masmol
   END FUNCTION LHeatX_sc
-    
+
   FUNCTION LHeatX_ve(temp,xESP) RESULT(res)
     !! Compute latent heat of a given specie at given temperature (vector).
@@ -208,10 +225,10 @@
     TYPE(mm_esp), INTENT(in)                   :: xESP !! Specie properties.
     REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Latent heat of given specie at given temperatures (\(J.kg^{-1}\)).
-    REAL(kind=mm_wp) :: ftm 
+    REAL(kind=mm_wp) :: ftm
     INTEGER      :: i
     DO i=1,SIZE(temp)
-      ftm=MIN(1._mm_wp-temp(i)/xESP%tc,1.e-3_mm_wp)
+      ftm=MAX(1._mm_wp-temp(i)/xESP%tc,1.e-3_mm_wp)
       res(i) = mm_rgas*xESP%tc*(7.08_mm_wp*ftm**0.354_mm_wp+10.95_mm_wp*xESP%w*ftm**0.456_mm_wp) / &
-               xESP%masmol 
+        xESP%masmol
     ENDDO
   END FUNCTION LHeatX_ve
@@ -219,5 +236,5 @@
   FUNCTION sigX_sc(temp,xESP) RESULT(res)
     !! Get the surface tension between a given specie and the air (scalar).
-    !! 
+    !!
     !! The method computes the surface tension equation as given in \cite{reid1986} p. 637 (eq. 12-3.6).
     REAL(kind=mm_wp), INTENT(in) :: temp !! temperature (K).
@@ -229,7 +246,7 @@
     sig = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
     sig = xESP%pc**(2._mm_wp/3._mm_wp)*xESP%tc**(1._mm_wp/3._mm_wp)*sig*(1._mm_wp-tr)**(11._mm_wp/9._mm_wp)
-    res = sig*1e3_mm_wp ! dyn/cm2 -> N/m
+    res = sig*1e-3_mm_wp ! dyn/cm -> N/m
   END FUNCTION sigX_sc
-  
+
   FUNCTION sigX_ve(temp,xESP) RESULT(res)
     !! Get the surface tension between a given specie and the air (vector).
@@ -240,11 +257,11 @@
     REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: res     !! Surface tensions (\(N.m^{-1}\)).
     INTEGER      :: i
-    REAL(kind=mm_wp) :: tr,tbr,sig
+    REAL(kind=mm_wp) :: tr,tbr,sig0,sig
     tbr = xESP%tb/xESP%tc
-    sig = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
+    sig0 = 0.1196_mm_wp*(1._mm_wp+(tbr*dlog(xESP%pc/1.01325_mm_wp))/(1._mm_wp-tbr))-0.279_mm_wp
     DO i=1,SIZE(temp)
       tr     = MIN(temp(i)/xESP%tc,0.99_mm_wp)
-      sig    = xESP%pc**(2._mm_wp/3._mm_wp)*xESP%tc**(1._mm_wp/3._mm_wp)*sig*(1._mm_wp-tr)**(11._mm_wp/9._mm_wp)
-      res(i) = sig*1e3_mm_wp ! dyn/cm2 -> N/m
+      sig    = xESP%pc**(2._mm_wp/3._mm_wp)*xESP%tc**(1._mm_wp/3._mm_wp)*sig0*(1._mm_wp-tr)**(11._mm_wp/9._mm_wp)
+      res(i) = sig*1e-3_mm_wp ! dyn/cm -> N/m
     ENDDO
   END FUNCTION sigX_ve
@@ -252,11 +269,6 @@
   FUNCTION psatX_sc(temp,xESP) RESULT(res)
     !! Get saturation vapor pressure for a given specie at given temperature (scalar).
-    !! 
+    !!
     !! The method computes the saturation vapor pressure equation given in \cite{reid1986} p. 657 (eq. 1).
-    !!
-    !! @warning
-    !! This subroutine accounts for a specific Titan feature: 
-    !! If __xESP__ corresponds to \(CH_{4}\), the saturation vapor presure is multiplied by 0.85 
-    !! to take into account its dissolution in \(N_{2}\).
     REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
     TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
@@ -266,14 +278,9 @@
     IF (x > 0._mm_wp) THEN
       qsat = (1._mm_wp-x)**(-1) * &
-      (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**3 + xESP%d_sat*x**6)
+        (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**2.5_mm_wp + xESP%d_sat*x**5_mm_wp)
     ELSE
-      qsat = XESP%a_sat*x/abs(1._mm_wp-x)     ! approx for  t > tc
-    ENDIF
-    !  Special case : ch4 : x0.85 (dissolution in N2)
-    IF (xESP%name == "ch4") THEN
-      res = xESP%pc*dexp(qsat)*0.85_mm_wp
-    ELSE
-      res = xESP%pc*dexp(qsat)
-    ENDIF
+      qsat = XESP%a_sat*x/abs(1._mm_wp-x) ! approx for t > tc
+    ENDIF
+    res = xESP%pc*exp(qsat)
     ! now convert bar to Pa
     res = res * 1e5_mm_wp
@@ -282,5 +289,5 @@
   FUNCTION psatX_ve(temp,xESP) RESULT(res)
     !! Get saturation vapor pressure for a given specie at given temperature (vector).
-    !! 
+    !!
     !! See [[mm_methods(module):psatX_sc(function)]].
     REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
@@ -292,29 +299,41 @@
       x = 1._mm_wp-temp(i)/xESP%tc
       IF (x > 0._mm_wp) THEN
-        qsat = (1._mm_wp-x)**(-1) * & 
-        (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**3 + xESP%d_sat*x**6)
+        qsat = (1._mm_wp-x)**(-1) * &
+          (xESP%a_sat*x + xESP%b_sat*x**1.5_mm_wp + xESP%c_sat*x**2.5_mm_wp + xESP%d_sat*x**5_mm_wp)
       ELSE
-        qsat = XESP%a_sat*x/abs(1._mm_wp-x)     ! approx for  t > tc
+        qsat = XESP%a_sat*x/abs(1._mm_wp-x) ! approx for t > tc
       ENDIF
-      res(i) = xESP%pc*dexp(qsat)
-      !  Peculiar case : ch4 : x0.85 (dissolution in N2)
-      IF (xESP%name == "ch4") res(i) = res(i)* 0.85_mm_wp 
+      res(i) = xESP%pc*exp(qsat)
     ENDDO
-    res = res * 1e5_mm_wp  ! bar -> Pa
+    ! now convert bar to Pa
+    res = res * 1e5_mm_wp
   END FUNCTION psatX_ve
 
   FUNCTION qsatX_sc(temp,pres,xESP) RESULT(res)
     !! Get the mass mixing ratio of a given specie at saturation (scalar).
+    !!
+    !! @warning
+    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
+    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
     REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
     REAL(kind=mm_wp), INTENT(in) :: pres !! Pressure level (Pa).
     TYPE(mm_esp), INTENT(in)    :: xESP  !! Specie properties.
     REAL(kind=mm_wp) :: res              !! Mass mixing ratio of the specie.
-    REAL(kind=mm_wp) :: x,psat
+    REAL(kind=mm_wp) :: psat
     psat = mm_psatX(temp,xESP)
     res = (psat / pres) * xESP%fmol2fmas
+    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
+    IF (xESP%name == "CH4") THEN
+      res = res * 0.85_mm_wp
+      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
+    ENDIF
   END FUNCTION qsatX_sc
 
   FUNCTION qsatX_ve(temp,pres,xESP) RESULT(res)
     !! Get the mass mixing ratio of a given specie at saturation (vector).
+    !!
+    !! @warning
+    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
+    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
     REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
     REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres !! Pressure levels (Pa).
@@ -324,5 +343,69 @@
     psat = mm_psatX(temp,xESP)
     res = (psat / pres) * xESP%fmol2fmas
+    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
+    IF (xESP%name == "CH4") THEN
+      res = res * 0.85_mm_wp
+      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
+    ENDIF
   END FUNCTION qsatX_ve
+
+  FUNCTION ysatX_sc(temp,pres,xESP) RESULT(res)
+    !! Get the molar mixing ratio of a given specie at saturation (scalar).
+    !!
+    !! @warning
+    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
+    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
+    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
+    REAL(kind=mm_wp), INTENT(in) :: pres !! Pressure level (Pa).
+    TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
+    REAL(kind=mm_wp)             :: res  !! Molar mixing ratio of the specie.
+
+    ! Fray and Schmidt (2009)
+    IF(xESP%name == "C2H2") THEN
+      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
+    
+    ELSE IF(xESP%name == "C2H6") THEN
+      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
+    
+    ELSE IF(xESP%name == "HCN") THEN
+      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
+    
+    ELSE IF (xESP%name == "CH4") THEN
+      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
+      res = res * 0.85_mm_wp
+      !IF (res < 0.014) THEN
+      !  res = 0.014
+      !ENDIF
+    ENDIF
+  END FUNCTION ysatX_sc
+
+  FUNCTION ysatX_ve(temp,pres,xESP) RESULT(res)
+    !! Get the molar mixing ratio of a given specie at saturation (vector).
+    !!
+    !! @warning
+    !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
+    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
+    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
+    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres !! Pressure levels (Pa).
+    TYPE(mm_esp), INTENT(in)                   :: xESP !! Specie properties.
+    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Molar mixing ratios of the specie.
+
+    ! Fray and Schmidt (2009)
+    IF(xESP%name == "C2H2") THEN
+      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
+    
+    ELSE IF(xESP%name == "C2H6") THEN
+      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
+    
+    ELSE IF(xESP%name == "HCN") THEN
+      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
+    
+    ! Peculiar case : CH4 : x0.85 (dissolution in N2)
+    ELSE IF (xESP%name == "CH4") THEN
+      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
+      res = res * 0.85_mm_wp
+      !WHERE (res(:) < 0.014) res(:) = 0.014
+    ENDIF
+  END FUNCTION ysatX_ve
 
   ELEMENTAL FUNCTION mm_get_kco(t) RESULT(res)
