diff --git a/source/model/numerics_mod.f90 b/source/model/numerics_mod.f90
index acdd4c3c1bbc24e297c19fd756dd36d1d68506a7..23ea5ae54704a4274e8e42ca6e0cd8412591b51a 100644
--- a/source/model/numerics_mod.f90
+++ b/source/model/numerics_mod.f90
@@ -350,4 +350,147 @@ y(1:N) = x(1:N)
 deallocate(x)
 END SUBROUTINE SMOOTHER_CONSERV
 
+!************************************************************************
+!*                                                                      *
+!*    Program to calculate the first kind modified Bessel function of   *
+!*    integer order N, for any REAL X, using the function BESSI(N,X).   *
+!*                                                                      *
+!* -------------------------------------------------------------------- *
+!*                                                                      *
+!*    SAMPLE RUN:                                                       *
+!*                                                                      *
+!*    (Calculate Bessel function for N=2, X=0.75).                      *
+!*                                                                      *
+!*    Bessel function of order  2 for X =  0.7500:                      *
+!*                                                                      *
+!*         Y =  0.73666878E-01                                          *
+!*                                                                      *
+!* -------------------------------------------------------------------- *
+!*   Reference: From Numath Library By Tuan Dang Trong in Fortran 77.   *
+!*                                                                      *
+!*                               F90 Release 1.2 By J-P Moreau, Paris.  *
+!*                                        (www.jpmoreau.fr)             *
+!*                                                                      *
+!*   Version 1.1: corected value of P4 in BESSIO (P4=1.2067492 and not  *
+!*                1.2067429) Aug. 2011.                                 *
+!*   Version 2: all variables are declared.                             *
+!************************************************************************
+!PROGRAM TBESSI
+!
+!  IMPLICIT NONE
+!  REAL*8  BESSI, X, Y
+!  INTEGER N
+!
+!  N=2
+!  X=0.75D0
+!
+!  Y = BESSI(N,X)
+!
+!  write(*,10) N, X
+!  write(*,20) Y
+!
+!  stop
+!
+!10 format (/' Bessel Function of order ',I2,' for X=',F8.4,':')
+!20 format(/'      Y = ',E15.8/)
+!
+!END
+
+! ----------------------------------------------------------------------
+! Auxiliary Bessel functions for N=0, N=1
+      FUNCTION BESSI0(X)
+      IMPLICIT NONE
+      REAL *8 X,BESSI0,Y,P1,P2,P3,P4,P5,P6,P7,  &
+      & Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
+      DATA P1,P2,P3,P4,P5,P6,P7/1.D0,3.5156229D0,3.0899424D0,1.2067492D0,  &
+      & 0.2659732D0,0.360768D-1,0.45813D-2/
+      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, &
+      & 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,  &
+      & 0.2635537D-1,-0.1647633D-1,0.392377D-2/
+      IF(ABS(X).LT.3.75D0) THEN
+      Y=(X/3.75D0)**2
+      BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
+      ELSE
+      AX=ABS(X)
+      Y=3.75D0/AX
+      BX=EXP(AX)/SQRT(AX)
+      AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
+      BESSI0=AX*BX
+      ENDIF
+      RETURN
+      END
+! ----------------------------------------------------------------------
+      FUNCTION BESSI1(X)
+      IMPLICIT NONE
+      REAL *8 X,BESSI1,Y,P1,P2,P3,P4,P5,P6,P7,  &
+      & Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX
+      DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,  &
+      & 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
+      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1, &
+      & -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, &
+      & -0.2895312D-1,0.1787654D-1,-0.420059D-2/
+      IF(ABS(X).LT.3.75D0) THEN
+      Y=(X/3.75D0)**2
+      BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
+      ELSE
+      AX=ABS(X)
+      Y=3.75D0/AX
+      BX=EXP(AX)/SQRT(AX)
+      AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
+      BESSI1=AX*BX
+      ENDIF
+      RETURN
+      END
+
+! ----------------------------------------------------------------------
+
+! ----------------------------------------------------------------------
+      FUNCTION BESSI(N,X)
+!
+!     This subroutine calculates the first kind modified Bessel function
+!     of integer order N, for any REAL X. We use here the classical
+!     recursion formula, when X > N. For X < N, the Miller's algorithm
+!     is used to avoid overflows. 
+!     REFERENCE:
+!     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
+!     MATHEMATICAL TABLES, VOL.5, 1962.
+!
+      IMPLICIT NONE
+      INTEGER, PARAMETER :: IACC = 40
+      REAL*8, PARAMETER :: BIGNO = 1.D10, BIGNI = 1.D-10
+      INTEGER N, M, J
+      REAL *8 X,BESSI,TOX,BIM,BI,BIP !,BESSI0,BESSI1
+      IF (N.EQ.0) THEN
+      BESSI = BESSI0(X)
+      RETURN
+      ENDIF
+      IF (N.EQ.1) THEN
+      BESSI = BESSI1(X)
+      RETURN
+      ENDIF
+      IF(X.EQ.0.D0) THEN
+      BESSI=0.D0
+      RETURN
+      ENDIF
+      TOX = 2.D0/X
+      BIP = 0.D0
+      BI  = 1.D0
+      BESSI = 0.D0
+      M = 2*((N+INT(SQRT(FLOAT(IACC*N)))))
+      DO 12 J = M,1,-1
+      BIM = BIP+DFLOAT(J)*TOX*BI
+      BIP = BI
+      BI  = BIM
+      IF (ABS(BI).GT.BIGNO) THEN
+      BI  = BI*BIGNI
+      BIP = BIP*BIGNI
+      BESSI = BESSI*BIGNI
+      ENDIF
+      IF (J.EQ.N) BESSI = BIP
+12    CONTINUE
+      BESSI = BESSI*BESSI0(X)/BI
+      RETURN
+      END
+
+
 END MODULE NUMERICS
diff --git a/source/model/phys_func.f90 b/source/model/phys_func.f90
index 317e5ab729a078e3dc09d620da8ca489b58aaa99..d5a6d3dd24535d3b65fcd10027d2aad10c5b34b0 100644
--- a/source/model/phys_func.f90
+++ b/source/model/phys_func.f90
@@ -1985,5 +1985,25 @@ END FUNCTION FP_THETA
   if (firstcall) firstcall = .false.
   END SUBROUTINE DIFFMIN_HS
 
+  !> Function HORIZCONC returns a value of dissolved concentration in the mixed layer
+  !! at distance r from center of circular lake; according to 
+  !! (DelSontro et al., Ecosystems, 2018) extended by inclusion of linear sink term
+  FUNCTION HORIZCONC(Caver,r,rL,kdiff,kexch,ksink,hML) result(Cr)
+  use NUMERICS, only : BESSI !Bessel function
+  implicit none
+  !Input/output variables
+  real(kind=ireals), intent(in) :: Caver !> horizontally averaged concentration in the ML
+  real(kind=ireals), intent(in) :: r     !> distance form lake center, m
+  real(kind=ireals), intent(in) :: rL    !> lake radius, m
+  real(kind=ireals), intent(in) :: kdiff !> effective diffusion coefficient in horizontal, m**2/s
+  real(kind=ireals), intent(in) :: kexch !> gas exchange koefficient with the atmosphere (piston velocity), m/s
+  real(kind=ireals), intent(in) :: ksink !> the rate of concentration sink in the ML, s**(-1)
+  real(kind=ireals), intent(in) :: hML   !> mixed-layer (ML) depth
+  real(kind=ireals)             :: Cr    !> concentration at the distance r from lake center
+  !Local variables
+  real(kind=ireals) :: x
+  x = sqrt((kexch/hML+ksink)/kdiff)
+  Cr = 0.5*Caver*rL*BESSI(0_4,x*r)/BESSI(1_4,x*rL)
+  END FUNCTION HORIZCONC
 
   END MODULE PHYS_FUNC