Skip to content
Snippets Groups Projects
Commit 9b67413a authored by Victor Stepanenko's avatar Victor Stepanenko
Browse files

Reconstruction of horizonal concentration distribution in a mixed layer of a circular lake is added

parent d60cf50a
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
source/model/numerics_mod.f90
+ 143
− 0
View file @ 9b67413a
......@@ -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
source/model/phys_func.f90
+ 20
− 0
View file @ 9b67413a
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment