Skip to content
Snippets Groups Projects
Commit 3d1309a5 authored by Ramil Ahtamyanov's avatar Ramil Ahtamyanov
Browse files

rename ocevermix module to obl_inmom. Separating the obl_legacy module into a separate file

parent 7b7685ea
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
obl_inmom.f90 0 → 100644
+ 161
− 0
View file @ 3d1309a5
module obl_inmom
implicit none
integer, parameter :: border_shift = 1
private :: border_shift
contains
subroutine init_vermix(richnum_mode, kh_km_mode, uu, vv, lu, dx, dy, dxh, dyh, hhu, hhv, hhq, zw, g, rh0, &
rit, den, rlh, taux, tauy, aice0, tt, ss, anzt, anzu, anumaxt, anumaxu, anubgrt, anubgru)
use obl_legacy
use obl_pph
use obl_pph_dyn
integer, intent(in) :: richnum_mode, kh_km_mode
real, intent(inout) :: anzt(:,:,:) ! kh
real, intent(inout) :: anzu(:,:,:) ! km
real, intent(inout) :: rit(:,:,:)
real, intent(inout) :: den(:,:,:)
real, intent(in) :: uu(:,:,:), vv(:,:,:)
real, intent(in) :: tt(:,:,:)
real, intent(in) :: ss(:,:,:)
real, intent(in) :: taux(:,:), tauy(:,:)
real, intent(in) :: lu(:,:)
real, intent(in) :: dx(:,:), dy(:,:), dxh(:,:), dyh(:,:)
real, intent(in) :: hhu(:,:), hhv(:,:), hhq(:,:)
real, intent(in) :: aice0(:,:)
real, intent(in) :: rlh(:,:)
real, intent(in) :: zw(:)
real, intent(in) :: g, rh0
real, intent(in) :: anumaxt ! kh_0
real, intent(in) :: anumaxu ! km_0
real, intent(in) :: anubgrt ! kh_b0
real, intent(in) :: anubgru ! km_b0
! Local variables
real, allocatable :: u2(:,:,:), v2(:,:,:)
real, allocatable :: s2(:,:,:)
real, allocatable :: n2(:,:,:)
real, allocatable :: kh_b(:)
real, allocatable :: km_b(:)
real, allocatable :: kh(:,:,:)
real, allocatable :: km(:,:,:)
real :: kh_str, km_str
real :: kh_undimdepth, km_undimdepth
real :: kh_unstable
real :: kh_b0, km_b0
real :: kh_0, km_0
integer :: nx, ny, nz
integer :: i, j
! Modules types
type(pphParamType) :: pphParams
type(gridDataType) :: inmomGrid
allocate(u2(size(uu, 1), size(uu, 2), size(uu, 3)))
allocate(v2(size(vv, 1), size(vv, 2), size(vv, 3)))
allocate(s2(size(u2, 1), size(u2, 2), size(u2, 3) - 1))
allocate(n2(size(den, 1), size(den, 2), size(den, 3) - 1))
allocate(kh(size(anzt, 1), size(anzt, 2), size(anzt, 3)))
allocate(km(size(anzu, 1), size(anzu, 2), size(anzu, 3)))
allocate(kh_b(size(kh, 3)))
allocate(km_b(size(km, 3)))
kh = anzt
km = anzu
kh_0 = anumaxt
km_0 = anumaxu
kh_b0 = anubgrt
km_b0 = anubgru
nx = size(rit, 1)
ny = size(rit, 2)
nz = size(rit, 3)
! print *, "hello 1"
den = legacy_denp(tt, ss + 35.0, 0.0)
if (richnum_mode == 1) then
! print *, "hello rit calc"
call legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
call legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
call legacy_s2(u2, v2, lu, s2)
call legacy_n2(den, hhq, zw, g, lu, n2)
call legacy_rit(n2, s2, border_shift, lu, rit(:,:,2:size(rit, 3)))
call legacy_rit_top(rlh, taux, tauy, border_shift, lu, rit(:,:,1))
end if
! print *, "hello aft rit"
call sync_xy_border(rit)
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
if (kh_km_mode == 1) then
call legacy_str("kh", taux(i,j), tauy(i,j), rh0, kh_str)
call legacy_undimdepth("kh", kh_str, hhq(i,j), aice0(i,j), kh_undimdepth)
call legacy_kh_b(zw(:), hhq(i,j), kh_unstable, kh_undimdepth, kh_b0, kh_b)
call legacy_kh_unstable(tt(i,j,1), kh_unstable)
call legacy_kh(rit(i,j,:), kh_0, kh_b(:), kh_unstable, kh(i,j,:))
call legacy_str("km", taux(i,j), tauy(i,j), rh0, km_str)
call legacy_undimdepth("km", km_str, hhq(i,j), aice0(i,j), km_undimdepth)
call legacy_km_b(zw, km_unstable, km_undimdepth, km_b0, km_b)
call legacy_km(rit(i,j,:), km_0, km_b(:), km(i,j,:))
else if (kh_km_mode == 2) then
! call get_Kh(kh(i,j,:), rit(i,j,:), nz)
! call get_Km(km(i,j,:), rit(i,j,:), nz)
call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), pphParams, inmomGrid)
call get_eddy_viscosity(km(i,j,:), rit(i,j,:), pphParams, inmomGrid)
else if (kh_km_mode == 3) then
! u_dynH =
! mld =
call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
call get_eddy_viscosity(km(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
end if
end if
end do
end do
if (kh_km_mode == 2) then
kh = kh * 10000.0
km = km * 10000.0
end if
! print *, "Kh first:", kh(3,3,1:4)
! print *, "Km first:", km(3,3,1:4)
! print *, "Kh last:", kh(3,3,size(kh,3)-2:size(kh,3))
! print *, "Km last:", km(3,3,size(km,3)-2:size(km,3))
anzt = kh
anzu = km
end subroutine init_vermix
subroutine sync_xy_border(array)
implicit none
real, intent(inout) :: array(:,:,:)
integer :: nx, ny, nz, k
! Determine array dimensions
nx = size(array, 1)
ny = size(array, 2)
nz = size(array, 3)
! Update boundary points along x (first_x and end_x boundaries)
do k = 1, nz
array(1,2:ny-1,k) = array(2,2:ny-1,k) ! first_x boundary
array(nx,2:ny-1,k) = array(nx-1,2:ny-1,k) ! end_x boundary
end do
! Update boundary points along y (first_y and end_y boundaries)
do k = 1, nz
array(2:nx-1,1,k) = array(2:nx-1,2,k) ! first_y boundary
array(2:nx-1,ny,k) = array(2:nx-1,ny-1,k) ! end_y boundary
end do
! Update corner points
do k = 1, nz
array(1,1,k) = (array(2,1,k) + array(1,2,k)) / 2.0
array(1,ny,k) = (array(2,ny,k) + array(1,ny-1,k)) / 2.0
array(nx,1,k) = (array(nx-1,1,k) + array(nx,2,k)) / 2.0
array(nx,ny,k) = (array(nx-1,ny,k) + array(nx,ny-1,k)) / 2.0
end do
end subroutine sync_xy_border
end module obl_inmom
\ No newline at end of file
obl_legacy.f90 0 → 100644
+ 465
− 0
View file @ 3d1309a5
module obl_legacy
implicit none
!!! OCEVERMIX LEGACY PARAMETERS
real, parameter :: lu_min = 0.5
! legacy_rit_top
real, parameter :: coriolis_rit_top_min = 2.5e-05
! legacy_rit
real, parameter :: rit_max = 1E+05
real, parameter :: s2_add = 0.0025
! legacy_kh_unstable
real, parameter :: kh_ttmx = 10.0
real, parameter :: kh_ttmn = 5.0
real, parameter :: kh_unstmx = 500.0
real, parameter :: kh_unstmn = 500.0 !50
! legacy_undimdepth & legacy_str
real, parameter :: hm_ice_max = 0.7
real, parameter :: hhq_add = 1.0e-30
real, parameter :: kh_strmin1 = 0.5
real, parameter :: kh_strmax1 = 2.0
real, parameter :: kh_strmin2 = 2.0
real, parameter :: kh_strmax2 = 10.0
real, parameter :: kh_hmmin1 = 300.0
real, parameter :: kh_hmmax1 = 3000.0
real, parameter :: kh_hmmin2 = 3000.0
real, parameter :: kh_hmmax2 = 10000.0
real, parameter :: kh_dimdepth = 1000.0
real, parameter :: km_strmin1 = 0.4
real, parameter :: km_strmax1 = 1.0
real, parameter :: km_strmin2 = 1.0
real, parameter :: km_strmax2 = 8.0
real, parameter :: km_hmmin1 = 300.0
real, parameter :: km_hmmax1 = 3000.0
real, parameter :: km_hmmin2 = 3000.0
real, parameter :: km_hmmax2 = 10000.0
real, parameter :: km_dimdepth = 1000.0 !1000.0
! legacy_kh_b
real, parameter :: kh_amnazmin = 1.0
real, parameter :: kh_amnazmax = 20.0
real, parameter :: kh_zzmin = 100000.0
real, parameter :: kh_zzmax = 150000.0
! legacy_km_b
real, parameter :: km_unstable = 500.0
!!! OCEVERMIX PARAMETERS VISIBILITY
public :: lu_min, km_unstable
private :: coriolis_rit_top_min, rit_max, s2_add, kh_ttmx, kh_ttmn, kh_unstmx, kh_unstmn
private :: hm_ice_max, hhq_add, kh_strmin1, kh_strmax1, kh_strmin2, kh_strmax2
private :: kh_hmmin1, kh_hmmax1, kh_hmmin2, kh_hmmax2, kh_dimdepth
private :: km_strmin1, km_strmax1, km_strmin2, km_strmax2
private :: km_hmmin1, km_hmmax1, km_hmmin2, km_hmmax2, km_dimdepth
private :: kh_amnazmin, kh_amnazmax, kh_zzmin, kh_zzmax
contains
subroutine legacy_rit(n2, s2, border_shift, lu, rit)
real, intent(in) :: n2(:,:,:)
real, intent(in) :: s2(:,:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: rit(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
! Parameters: rit_max, s2_add, lu_min
nx = size(rit, 1)
ny = size(rit, 2)
nz = size(rit, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
rit(i,j,k) = n2(i,j,k) / (s2(i,j,k) + s2_add)
rit(i,j,k) = min(rit(i,j,k), rit_max)
end if
end do
end do
end do
end subroutine legacy_rit
subroutine legacy_rit_top(rlh, tau_u, tau_v, border_shift, lu, rit_top)
real, intent(in) :: rlh(:,:)
real, intent(in) :: tau_u(:,:), tau_v(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: rit_top(:,:)
! Local variables
real :: rlt, coriolis, u_star ! Auxiliary variables
integer :: i, j ! Loop indices
integer :: nx, ny ! Array sizes
! Parameters: coriolis_rit_top_min, lu_min
nx = size(rit_top, 1)
ny = size(rit_top, 2)
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
rlt = (rlh(i,j) + rlh(i-1,j) + rlh(i,j-1) + rlh(i-1,j-1)) / 4.0
coriolis = max(abs(rlt), coriolis_rit_top_min) ! freezing in 10N,S
u_star = sqrt(sqrt(tau_u(i,j)**2 + tau_v(i,j)**2))
rit_top(i,j) = 0.5 * u_star / coriolis
end if
end do
end do
end subroutine legacy_rit_top
subroutine legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
real, intent(in) :: uu(:,:,:)
real, intent(in) :: dy(:,:)
real, intent(in) :: dyh(:,:)
real, intent(in) :: hhu(:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: u2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(uu, 1)
ny = size(uu, 2)
nz = size(uu, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
u2(i,j,k) = lu(i,j) * 0.5 / dy(i,j) * &
(uu(i-1,j,k) * hhu(i-1,j) * dyh(i-1,j) + &
uu(i,j,k) * hhu(i,j) * dyh(i,j)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_u2
subroutine legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
real, intent(in) :: vv(:,:,:)
real, intent(in) :: dx(:,:)
real, intent(in) :: dxh(:,:)
real, intent(in) :: hhv(:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
integer, intent(in) :: border_shift
real, intent(inout) :: v2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(vv, 1)
ny = size(vv, 2)
nz = size(vv, 3)
do k = 1, nz
do j = 1 + border_shift, ny - border_shift
do i = 1 + border_shift, nx - border_shift
if (lu(i,j) > lu_min) then
v2(i,j,k) = lu(i,j) * 0.5 / dx(i,j) * &
(vv(i,j-1,k) * hhv(i,j-1) * dxh(i,j-1) + &
vv(i,j,k) * hhv(i,j) * dxh(i,j)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_v2
subroutine legacy_s2(u2, v2, lu, s2)
real, intent(in) :: u2(:,:,:)
real, intent(in) :: v2(:,:,:)
real, intent(in) :: lu(:,:)
real, intent(inout) :: s2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(s2, 1)
ny = size(s2, 2)
nz = size(s2, 3)
do k = 1, nz
do j = 1, ny
do i = 1, nx
if (lu(i,j) > lu_min) then
s2(i,j,k) = (u2(i,j,k+1) - u2(i,j,k))**2 + (v2(i,j,k+1) - v2(i,j,k))**2
end if
end do
end do
end do
end subroutine legacy_s2
subroutine legacy_n2(den, hhq, zw, g, lu, n2)
real, intent(in) :: den(:,:,:)
real, intent(in) :: hhq(:,:)
real, intent(in) :: lu(:,:)
real, intent(in) :: zw(:)
real, intent(in) :: g
real, intent(inout) :: n2(:,:,:)
! Local variables
integer :: i, j, k ! Loop indices
integer :: nx, ny, nz ! Array sizes
nx = size(n2, 1)
ny = size(n2, 2)
nz = size(n2, 3)
do k = 1, nz
do j = 1, ny
do i = 1, nx
if (lu(i,j) > lu_min) then
n2(i,j,k) = g * (den(i,j,k+1) - den(i,j,k)) / 2.0 * &
(zw(k+2) - zw(k)) / hhq(i,j)
end if
end do
end do
end do
end subroutine legacy_n2
subroutine legacy_km_b(zw, unstable, undimdepth, km_b0, km_b)
real, intent(in) :: zw(:)
real, intent(in) :: unstable
real, intent(in) :: undimdepth
real, intent(in) :: km_b0
real, intent(inout) :: km_b(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(km_b) - 1
do k = 2, nz
! if (zw(k) <= undimdepth) then
! km_b(k) = unstable
! else
! km_b(k) = km_b0
! end if
km_b(k) = km_b0
end do
end subroutine legacy_km_b
subroutine legacy_km(rit, km_0, km_b, km)
real, intent(in) :: rit(:)
real, intent(in) :: km_0
real, intent(in) :: km_b(:)
real, intent(inout) :: km(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(rit)
do k = 2, nz
if (rit(k) > 0.0) then
km(k) = km_0 / (1.0 + 5.0 * rit(k))**2 + km_b(k)
else
km(k) = km_0 + km_b(k)
end if
end do
end subroutine legacy_km
subroutine legacy_kh(rit, kh_0, kh_b, unstable, kh)
real, intent(in) :: rit(:)
real, intent(in) :: kh_0
real, intent(in) :: unstable
real, intent(in) :: kh_b(:)
real, intent(inout) :: kh(:)
! Local variables
integer :: k ! Loop indices
integer :: nz ! Array sizes
nz = size(rit)
do k = 2, nz
if (rit(k) >= 0.0) then
kh(k) = kh_0 / (1.0 + 5.0 * rit(k)) + kh_b(k)
else
kh(k) = unstable
end if
end do
! For unstable stratification
! do k = 2, nz
! if (rit(k) < 0.0) then
! Mixing is decreased in vicinity of equator in unstable situation
! kh(k) = unstable * (1.0 - 0.7 * (rn / rm) ** 10)
! kh(k) = unstable
! kh(k + 1) = unstable
! ! Convective diffusion as dZ * dZ * Vajsjala-Brenta frequency
! kh(k) = az0 * hhq * hzt(k) * &
! sqrt(abs(den(k) - den(k - 1)) * &
! g * hhq * hzt(k)) + &
! 0.3 * hhq**2 * hzt(k-1) * &
! hzt(k) / (12.0 * 3600.0)
! end if
! end do
end subroutine legacy_kh
subroutine legacy_kh_b(zw, hhq, unstable, undimdepth, kh_b0, kh_b)
real, intent(in) :: zw(:)
real, intent(in) :: hhq
real, intent(in) :: unstable
real, intent(in) :: undimdepth
real, intent(in) :: kh_b0
real, intent(inout) :: kh_b(:)
! Local variables
real :: amnaz, zz
integer :: k ! Loop indices
integer :: nz ! Array sizes
! Parameters: kh_amnazmin, kh_amnazmax, kh_zzmin, kh_zzmax
nz = size(kh_b) - 1
do k = 2, nz
if (zw(k) <= undimdepth) then
kh_b = unstable ! Addition mixing in upper layer for T & S
else
! BACKGROUND COEFFICIENT DEPENDS ON DEPTH
zz = zw(k) * hhq
if (zz <= kh_zzmin) then
amnaz = kh_amnazmin
else if (zz >= kh_zzmax) then
amnaz = kh_amnazmax
else
amnaz = kh_amnazmin + (kh_amnazmax - kh_amnazmin) * &
(zz - kh_zzmin) / (kh_zzmax - kh_zzmin)
end if
kh_b(k) = kh_b0 * amnaz
end if
end do
end subroutine legacy_kh_b
subroutine legacy_undimdepth(mode, str, hhq, aice0, undimdepth)
character(len=2) :: mode
real, intent(in) :: str
real, intent(in) :: hhq
real, intent(in) :: aice0
real, intent(inout) :: undimdepth
! Local variables
real :: strmin1, strmax1, strmin2, strmax2
real :: hmmin1, hmmax1, hmmin2, hmmax2
real :: dimdepth
real :: hm
! Parameters: kh_strmin1, kh_strmax1, kh_strmin2, kh_strmax2
! kh_hmmin1, kh_hmmax1, kh_hmmin2, kh_hmmax2
! kh_dimdepth
! km_strmin1, km_strmax1, km_strmin2, km_strmax2
! km_hmmin1, km_hmmax1, km_hmmin2, km_hmmax2
! km_dimdepth
if (mode == "kh") then
strmin1 = kh_strmin1
strmax1 = kh_strmax1
strmin2 = kh_strmin2
strmax2 = kh_strmax2
hmmin1 = kh_hmmin1
hmmax1 = kh_hmmax1
hmmin2 = kh_hmmin2
hmmax2 = kh_hmmax2
dimdepth = kh_dimdepth
else if (mode == "km") then
strmin1 = km_strmin1
strmax1 = km_strmax1
strmin2 = km_strmin2
strmax2 = km_strmax2
hmmin1 = km_hmmin1
hmmax1 = km_hmmax1
hmmin2 = km_hmmin2
hmmax2 = km_hmmax2
dimdepth = km_dimdepth
else
stop "legacy_undimdepth: ERROR: mode must be 'kh' or 'km'"
end if
if (str <= strmax1) then
hm = hmmin1 + (hmmax1 - hmmin1) * (str - strmin1) / (strmax1 - strmin1)
else
hm = hmmin2 + (hmmax2 - hmmin2) * (str - strmin2) / (strmax2 - strmin2)
end if
if (aice0 < hm_ice_max) then
hm = dimdepth
end if
undimdepth = amax1(dimdepth, hm) / (hhq + hhq_add)
end subroutine legacy_undimdepth
subroutine legacy_str(mode, taux, tauy, rh0, str)
character(len=2) :: mode
real, intent(in) :: taux, tauy
real, intent(in) :: rh0
real, intent(inout) :: str
! Local variables
real :: strmin, strmax
! Parameters: kh_strmin1, kh_strmax2, km_strmin1, km_strmax2
if (mode == "kh") then
strmin = kh_strmin1
strmax = kh_strmax2
else if (mode == "km") then
strmin = km_strmin1
strmax = km_strmax2
else
stop "legacy_str: ERROR: mode must be 'kh' or 'km'"
end if
str = rh0 * sqrt(taux**2 + tauy**2)
str = amax1(str, strmin) ! TODO: check is amax1, amin1 is not outdated
str = amin1(str, strmax)
end subroutine legacy_str
subroutine legacy_kh_unstable(tt, unstable)
real, intent(in) :: tt
real, intent(inout) :: unstable
! Parameters: kh_ttmx, kh_ttmn, kh_unstmx, kh_unstmn
if (tt > kh_ttmx) then
unstable = kh_unstmx
else if (tt < kh_ttmn) then
unstable = kh_unstmn
else
unstable = kh_unstmn + (kh_unstmx - kh_unstmn) * &
(tt - kh_ttmn) / (kh_ttmx - kh_ttmn)
end if
end subroutine legacy_kh_unstable
!> Computes the sea water density deviation from 1.02 [g/cm³].
!!
!! This function calculates the potential density deviation
!! of seawater as a function of potential temperature (**t** in °C),
!! salinity (**s** in parts per thousand), and pressure (**p** in MPa).
!! real, parameter :: prestompa = 1.0e-7 ! conversion factor for pressure [1MPa=1E-7din/cm**2] <init: Inc/0DENP.INC; using in: occont.f: ocforcing>
!!
!! **References:**
!! The calculation based on Bryden et al. (1999): "A new approximation of
!! the equation of state for seawater, suitable for numerical ocean models."
!! https://doi.org/10.1029/1998JC900059
!!
!! **Valid ranges:**
!! - Temperature: -2 < t < 40 °C
!! - Salinity: 0 < s < 42 ppt
!! - Pressure: 0 < p < 100 MPa
!!
!! **Input Variables:**
!! @param[in] t Potential temperature in degrees Celsius.
!! @param[in] s Salinity in parts per thousand (ppt).
!! @param[in] p Pressure in megapascals (MPa).
!! **Return Value:**
!! @return Potential density deviation from 1.02 [g/cm³].
elemental real function legacy_denp(t, s, p) result(denp)
real, intent(in) :: t, s, p
denp = -2.00920601e-02 + &
( 5.07043e-04 * p - 5.43283e-07 * p * p) + &
( 5.10768e-05 - 3.69119e-06 * p + 6.54837e-09 * p * p) * t + &
( 8.05999e-04 - 9.34012e-07 * p + 1.38777e-09 * p * p) * s + &
(-7.40849e-06 + 5.33243e-08 * p - 1.01563e-10 * p * p) * t * t + &
(-3.01036e-06 + 1.75145e-08 * p - 2.34892e-11 * p * p) * t * s + &
( 3.32267e-08 - 3.25887e-10 * p + 4.98612e-13 * p * p) * t * t * t + &
( 3.21931e-08 - 1.65849e-10 * p + 2.17612e-13 * p * p) * t * t * s
end function legacy_denp
end module obl_legacy
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment