Skip to content

GitLab

U
UVic-updates-opem
Commit 749d30fc authored by Markus Pahlow's avatar Markus Pahlow
Browse files
Options

make code compilable by gfortran

parent 653ae5a6
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 370 additions and 370 deletions
f90/cfo.f90
View file @ 749d30fc
......@@ -3,36 +3,36 @@ MODULE cfo
IMPLICIT NONE
PRIVATE
PUBLIC :: ft_ME, nprey
INTEGER, PARAMETER :: dp=KIND(0D0)
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15, 60)
INTEGER :: nprey=4
REAL(dp), PARAMETER :: cfds=86.4D3
REAL(dp), PARAMETER :: cfds=86.4E3_dp
INCLUDE 'stdunits.h'
TYPE, PUBLIC :: ocf
REAL(dp) :: ngr,& ! relative growth rate
Imax=1D0,& ! max. specific ingestion rate (1/d)
Rm=0.05D0,& ! specific maintenance respiration (1/d)
beta=0.2D0,& ! digestion (assimilation) coefficient
ca=0.1D0, cf=0.1D0,& ! cost of assimilation, foraging coefficients
Emax=1D0,& ! max. assimilation efficiency
gbio, geZ, gamma1, kzoo,&
Af,& ! foraging activity
At,& ! total activity
IC, IN, IP,& ! C, N, P ingestion
E,& ! assimilation efficiency
g,& ! net growth rate
POC, PON, POP,& ! zooplankton biomass concentration
RC, RN, RP,& ! respiration
XC,& ! C excretion
XN,& ! N excretion
XP,& ! P excretion
QN=0.15D0,& ! N:C ratio (mol/mol)
QP=0.0095D0,& ! P:C ratio (mol/mol)
Pth,& ! feeding threshold
morq, mort ! quadratic mortality (m3/mmolC/d)
REAL(dp) :: ngr,& ! relative growth rate
Imax=1E0_dp,& ! max. specific ingestion rate (1/d)
Rm=0.05E0_dp,& ! specific maintenance respiration (1/d)
beta=0.2E0_dp,& ! digestion (assimilation) coefficient
ca=0.1E0_dp, cf=0.1E0_dp,& ! cost of assimilation, foraging coefficients
Emax=1E0_dp,& ! max. assimilation efficiency
gbio, geZ, gamma1, kzoo,& ! UVic parameters
Af,& ! foraging activity
At,& ! total activity
IC, IN, IP,& ! C, N, P ingestion
E,& ! assimilation efficiency
g,& ! net growth rate
POC, PON, POP,& ! zooplankton biomass concentration
RC, RN, RP,& ! respiration
XC,& ! C excretion
XN,& ! N excretion
XP,& ! P excretion
QN=0.15E0_dp,& ! N:C ratio (mol/mol)
QP=0.0095E0_dp,& ! P:C ratio (mol/mol)
Pth,& ! feeding threshold
morq, mort ! quadratic mortality (m3/mmolC/d)
REAL(dp), DIMENSION(:), ALLOCATABLE :: phi,& ! prey capture coefficients (m3/molC)
pref,&
pI ! partial ingestion of each prey
pI ! partial ingestion of each prey
PROCEDURE(REAL(dp)), POINTER, NOPASS :: fT
PROCEDURE(foraging), POINTER :: forage
PROCEDURE(loss), POINTER :: regex
......@@ -42,16 +42,16 @@ MODULE cfo
ABSTRACT INTERFACE
SUBROUTINE foraging (zoo, preyC, preyN, preyP, bct, bctz)
IMPORT :: ocf
IMPORT :: ocf, dp
IMPLICIT NONE
CLASS(ocf), INTENT(INOUT) :: zoo
REAL(KIND(0D0)), INTENT(IN) :: preyC(:), preyN(:), preyP(:), bct, bctz
REAL(dp), INTENT(IN) :: preyC(:), preyN(:), preyP(:), bct, bctz
END SUBROUTINE foraging
SUBROUTINE loss (zoo, Rm, fQ)
IMPORT :: ocf
IMPORT :: ocf, dp
IMPLICIT NONE
CLASS(ocf), INTENT(INOUT) :: zoo
REAL(KIND(0D0)), INTENT(IN) :: Rm, fQ
REAL(dp), INTENT(IN) :: Rm, fQ
END SUBROUTINE loss
END INTERFACE
......@@ -79,13 +79,13 @@ CONTAINS
ca = zoo%ca
cf = zoo%cf
Emax = zoo%Emax
phi = 0D0
phi = 0E0_dp
Imax = zoo%Imax
Rm = zoo%Rm
QN = zoo%QN
QP = zoo%QP
morq0 = zoo%morq
morq = -1D0
morq = -1E0_dp
ftz = 'bctz'
regex = 'respire'
forage = ''
......@@ -94,15 +94,15 @@ CONTAINS
zoo%ca = ca
zoo%cf = cf
zoo%Emax = Emax
zoo%phi = phi*1D-3 ! convert m3/molC to m3/mmolC
zoo%phi = phi*1E-3_dp ! convert m3/molC to m3/mmolC
zoo%Imax = Imax/cfds ! convert 1/d to 1/s
zoo%Rm = Rm/cfds ! convert 1/d to 1/s
zoo%QN = QN
zoo%QP = QP
zoo%morq = MAX(morq, 0D0)/cfds
zoo%morq = MAX(morq, 0E0_dp)/cfds
zoo%At = zoo%Imax/zoo%beta & ! total activity (15)
*(-1D0 - lwm1(-(1D0 - zoo%cf/zoo%Emax/(1D0 - zoo%ca))/EXP(1D0 + zoo%beta)))
zoo%Pth = -LOG((1D0 - zoo%cf/zoo%Emax/(1D0 - zoo%ca))) ! /zoo%phi
*(-1E0_dp - lwm1(-(1E0_dp - zoo%cf/zoo%Emax/(1E0_dp - zoo%ca))/EXP(1E0_dp + zoo%beta)))
zoo%Pth = -LOG((1E0_dp - zoo%cf/zoo%Emax/(1E0_dp - zoo%ca))) ! /zoo%phi
IF (LEN_TRIM(forage).EQ.0) forage = 'ocf'
100 CONTINUE
SELECT CASE (TRIM(ftz))
......@@ -117,11 +117,11 @@ CONTAINS
CASE ('ocf', 'cfo')
zoo%forage => forage_cfo
CASE ('UVic', '', 'H2', 'Holling 2', 'Holling II')
IF (morq.LT.0D0) zoo%morq = morq0*zoo%QN/cfds ! -> UVic default in m3/mmolN/s
IF (morq.LT.0E0_dp) zoo%morq = morq0*zoo%QN/cfds ! -> UVic default in m3/mmolN/s
zoo%gbio = zoo%gbio/cfds
zoo%forage => forage_UVic
CASE ('H3', 'Holling 3', 'Holling III')
IF (morq.LT.0D0) zoo%morq = morq0*zoo%QN/cfds ! -> UVic default in m3/mmolN/s
IF (morq.LT.0E0_dp) zoo%morq = morq0*zoo%QN/cfds ! -> UVic default in m3/mmolN/s
zoo%gbio = zoo%gbio/cfds
zoo%forage => forage_H3
CASE DEFAULT
......@@ -146,14 +146,14 @@ CONTAINS
REAL(dp) :: At, fT, I0, Rm, PC, PN, PP, fQ, rxm
fT = zoo%fT(bct, bctz)
Rm = zoo%Rm*fT
PC = MAX(SUM(zoo%phi*preyC), 1D-30) ! (effective food concentration)*phi
PC = MAX(SUM(zoo%phi*preyC), 1E-30_dp) ! (effective food concentration)*phi
PN = SUM(zoo%phi*preyN)
PP = SUM(zoo%phi*preyP)
IF ((PC.LE.zoo%Pth).OR.(fT.LE.0D0)) THEN ! food below feeding threshold
fQ = 1D0
zoo%Af = 0D0
IF ((PC.LE.zoo%Pth).OR.(fT.LE.0E0_dp)) THEN ! food below feeding threshold
fQ = 1E0_dp
zoo%Af = 0E0_dp
zoo%E = zoo%Emax
zoo%IC = 0D0
zoo%IC = 0E0_dp
ELSE
fQ = MIN(PN/zoo%QN, PP/zoo%QP, PC)/PC
! fQ accounts for different N:C ratios in predator and prey, after
......@@ -161,24 +161,24 @@ CONTAINS
! to Kiørboe (1989), so Af should not be affected either. Thus, food
! C:N should influence only R and ngr. The same concept is applied here
! also to C:P.
I0 = 1D0 - EXP(-PC) ! ingestion saturation (6)
I0 = 1E0_dp - EXP(-PC) ! ingestion saturation (6)
At = fT*zoo%At ! total activity (15)
! Af: foraging activity
zoo%Af = zoo%beta*At/(-1D0 - lwm1(-(1D0 - zoo%cf/(zoo%Emax*(1 - zoo%ca)*I0)) &
/EXP(1D0 + zoo%beta))); ! (12)
zoo%Af = zoo%beta*At/(-1E0_dp - lwm1(-(1E0_dp - zoo%cf/(zoo%Emax*(1 - zoo%ca)*I0)) &
/EXP(1E0_dp + zoo%beta)));! (12)
! E: assimilation efficiency
zoo%E = zoo%Emax*(1D0 - EXP(-zoo%beta*(At/zoo%Af - 1D0))) ! (8)
zoo%E = zoo%Emax*(1E0_dp - EXP(-zoo%beta*(At/zoo%Af - 1E0_dp)))! (8)
zoo%IC = zoo%POC*I0*zoo%Af ! C ingestion (7)
END IF
zoo%pI = zoo%IC*zoo%phi/PC ! ingestion of each prey type
zoo%IN = zoo%IC*PN/PC ! N intake
zoo%IP = zoo%IC*PP/PC ! P intake
CALL zoo%regex(Rm, fQ)
rxm = zoo%XC/MAX(zoo%RC, 1D-30) ! ratio of egestion : metabolic losses
zoo%RN = (zoo%IN - zoo%ngr*zoo%QN)/(1D0 + rxm) ! DIN excretion
zoo%RP = (zoo%IP - zoo%ngr*zoo%QP)/(1D0 + rxm) ! DIP excretion
zoo%XN = zoo%RN*rxm ! PON egestion
zoo%XP = zoo%RP*rxm ! POP egestion
rxm = zoo%XC/MAX(zoo%RC, 1E-30_dp) ! ratio of egestion : metabolic losses
zoo%RN = (zoo%IN - zoo%ngr*zoo%QN)/(1E0_dp + rxm) ! DIN excretion
zoo%RP = (zoo%IP - zoo%ngr*zoo%QP)/(1E0_dp + rxm) ! DIP excretion
zoo%XN = zoo%RN*rxm ! PON egestion
zoo%XP = zoo%RP*rxm ! POP egestion
END SUBROUTINE forage_cfo
SUBROUTINE forage_UVic (zoo, preyC, preyN, preyP, bct, bctz)
......@@ -197,14 +197,14 @@ CONTAINS
! geZ_eff is (growth efficiency)*zoo%IC to prevent division by 0
geZ_eff = MIN(zoo%geZ*zoo%IC, zoo%gamma1*zoo%IN/zoo%QN, zoo%gamma1*zoo%IP/zoo%QP)
! Rx fluxes go into dissolved inorganic pools
zoo%RC = MAX(zoo%gamma1*zoo%IC - geZ_eff, 0D0)
zoo%RN = MAX(zoo%gamma1*zoo%IN - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QN, 0D0)
zoo%RP = MAX(zoo%gamma1*zoo%IP - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QP, 0D0)
zoo%RC = MAX(zoo%gamma1*zoo%IC - geZ_eff, 0E0_dp)
zoo%RN = MAX(zoo%gamma1*zoo%IN - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QN, 0E0_dp)
zoo%RP = MAX(zoo%gamma1*zoo%IP - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QP, 0E0_dp)
zoo%ngr = zoo%gamma1*zoo%IC - zoo%RC
! Xx fluxes go to detritus
zoo%XC = (1D0 - zoo%gamma1)*zoo%IC
zoo%XN = (1D0 - zoo%gamma1)*zoo%IN
zoo%XP = (1D0 - zoo%gamma1)*zoo%IP
zoo%XC = (1E0_dp - zoo%gamma1)*zoo%IC
zoo%XN = (1E0_dp - zoo%gamma1)*zoo%IN
zoo%XP = (1E0_dp - zoo%gamma1)*zoo%IP
END SUBROUTINE forage_UVic
SUBROUTINE forage_H3 (zoo, preyC, preyN, preyP, bct, bctz)
......@@ -225,14 +225,14 @@ CONTAINS
! geZ_eff is (growth efficiency)*zoo%IC to prevent division by 0
geZ_eff = MIN(zoo%geZ*zoo%IC, zoo%gamma1*zoo%IN/zoo%QN, zoo%gamma1*zoo%IP/zoo%QP)
! Rx fluxes go into dissolved inorganic pools
zoo%RC = MAX(zoo%gamma1*zoo%IC - geZ_eff, 0D0)
zoo%RN = MAX(zoo%gamma1*zoo%IN - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QN, 0D0)
zoo%RP = MAX(zoo%gamma1*zoo%IP - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QP, 0D0)
zoo%RC = MAX(zoo%gamma1*zoo%IC - geZ_eff, 0E0_dp)
zoo%RN = MAX(zoo%gamma1*zoo%IN - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QN, 0E0_dp)
zoo%RP = MAX(zoo%gamma1*zoo%IP - (zoo%gamma1*zoo%IC - zoo%RC)*zoo%QP, 0E0_dp)
zoo%ngr = zoo%gamma1*zoo%IC - zoo%RC
! Xx fluxes go to detritus
zoo%XC = (1D0 - zoo%gamma1)*zoo%IC
zoo%XN = (1D0 - zoo%gamma1)*zoo%IN
zoo%XP = (1D0 - zoo%gamma1)*zoo%IP
zoo%XC = (1E0_dp - zoo%gamma1)*zoo%IC
zoo%XN = (1E0_dp - zoo%gamma1)*zoo%IN
zoo%XP = (1E0_dp - zoo%gamma1)*zoo%IP
END SUBROUTINE forage_H3
! excrete (respire) extra C in food
......@@ -242,9 +242,9 @@ CONTAINS
REAL(dp), INTENT(IN) :: Rm, fQ
REAL(dp) :: EI
EI = zoo%E*zoo%IC
zoo%ngr = (EI*(1D0 - zoo%ca) - zoo%poc*(zoo%cf*zoo%Af + Rm))*fQ
zoo%RC = EI - zoo%ngr ! respiration
zoo%XC = zoo%IC*(1D0 - zoo%E) ! egestion
zoo%ngr = (EI*(1E0_dp - zoo%ca) - zoo%poc*(zoo%cf*zoo%Af + Rm))*fQ
zoo%RC = EI - zoo%ngr ! respiration
zoo%XC = zoo%IC*(1E0_dp - zoo%E) ! egestion
END SUBROUTINE zxrex
! egest extra C in food
......@@ -264,7 +264,7 @@ CONTAINS
IMPLICIT NONE
REAL(dp), INTENT(IN) :: temperature
REAL(dp) :: fT
fT = (1D0 + (temperature - 15D0)*9.5D0)/260D0
fT = (1E0_dp + (temperature - 15E0_dp)*9.5E0_dp)/260E0_dp
END FUNCTION fT_ME
! temperature function ftz_bct: use bct also for zooplankton
......
f90/cmo.f90
View file @ 749d30fc
......@@ -3,13 +3,13 @@ MODULE cmo
USE lambert
IMPLICIT NONE
PRIVATE
INTEGER, PARAMETER :: dp=KIND(0D0)
REAL(dp), PARAMETER :: cfds=86400D0
INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15, 60)
REAL(dp), PARAMETER :: cfds=86400E0_dp
INCLUDE 'stdunits.h'
TYPE :: lost
REAL(dp) :: partC, partN, partP,& ! particulate losses
dissC=0D0, dissN=0D0, dissP=0D0 ! dissolved losses
REAL(dp) :: partC, partN, partP,& ! particulate losses
dissC=0E0_dp, dissN=0E0_dp, dissP=0E0_dp ! dissolved losses
END TYPE lost
! about the mortality functions:
......@@ -29,18 +29,18 @@ MODULE cmo
! units of alpha: 1 m2/E*mol/gChl = 0.4 m2/W*mol/gChl/d
TYPE :: ocm
CHARACTER(LEN=15) :: morfun='' ! mortality function
REAL(dp) :: A0=1D2,& ! nutrient affinity in m3/mol/d
alpha=0.4D0,& ! light affinity in m2/W*mol/gChl/d
REAL(dp) :: A0=1E2_dp,& ! nutrient affinity in m3/mol/d
alpha=0.4E0_dp,& ! light affinity in m2/W*mol/gChl/d
QsN,& ! structure-bound N quota
Q0N=0.04D0,& ! subsistence N quota
Q0P=0.001D0,& ! subsistence P quota
rdl=0.5D0,& ! day-length parameter
RC=0.1D0,& ! Chl maintenance respiration
V0=5D0,& ! maximum rate parameter
F0N=1.3D0,& ! maximum N2 fixation rate
zC=0.5D0,& ! cost of photosynthesis in mol/gChl
zN=0.6D0,& ! cost of biosynthesis in mol/mol
zF=2D0,& ! cost of N2 fixation
Q0N=0.04E0_dp,& ! subsistence N quota
Q0P=0.001E0_dp,& ! subsistence P quota
rdl=0.5E0_dp,& ! day-length parameter
RC=0.1E0_dp,& ! Chl maintenance respiration
V0=5E0_dp,& ! maximum rate parameter
F0N=1.3E0_dp,& ! maximum N2 fixation rate
zC=0.5E0_dp,& ! cost of photosynthesis in mol/gChl
zN=0.6E0_dp,& ! cost of biosynthesis in mol/mol
zF=2E0_dp,& ! cost of N2 fixation
QN, QNmax,& ! N quota (N:C ratio) in mol/mol
QP,& ! P quota (P:C ratio) in mol/mol
PARs, PARb,& ! irradiance
......@@ -56,15 +56,15 @@ MODULE cmo
FixN,& ! rate of N2 fixation
fV,& ! N allocation to nutrient acquisition
fN,& ! local N allocation to N acquisition
fF=0D0,& ! local N allocation to N2 fixation
fF=0E0_dp,& ! local N allocation to N2 fixation
fC,& ! N allocation to C fixation
SI,& ! degree of light saturation
tch,& ! local (chloroplast) Chl:C ratio
theta,& ! cellular Chl:C ratio
rctht,& ! rate of change of the Chl:C ratio
fcdtdq,& ! derivative fC/tht*d(tht)/d(QN)
mor=0D0, mort=0D0,& ! (quadratic) mortality, T-dependent mortality
lamor=0D0 ! Mayzaud-Poulet mortality coefficient in m3/mmol
mor=0E0_dp, mort=0E0_dp,& ! (quadratic) mortality, T-dependent mortality
lamor=0E0_dp ! Mayzaud-Poulet mortality coefficient in m3/mmol
REAL(dp), POINTER :: DIN,& ! DIN concentration
DIP ! DIP concentration
TYPE(lost) :: loss
......@@ -79,26 +79,26 @@ MODULE cmo
ABSTRACT INTERFACE
SUBROUTINE mort (phyt, phyc, bct)
IMPORT :: ocm
IMPORT :: ocm, dp
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(IN) :: phyc, bct
REAL(dp), INTENT(IN) :: phyc, bct
END SUBROUTINE mort
SUBROUTINE aisi (phyt, RC, phin, zN)
IMPORT :: ocm
IMPORT :: ocm, dp
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(INOUT) :: RC
REAL(KIND(0D0)), INTENT(IN) :: phin, zN
REAL(dp), INTENT(INOUT) :: RC
REAL(dp), INTENT(IN) :: phin, zN
END SUBROUTINE aisi
SUBROUTINE n2f (phy, felim, bct, zN) ! facultative N2 fixation
IMPORT :: ocm
IMPORT :: ocm, dp
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT) :: phy
REAL(KIND(0D0)), INTENT(IN) :: felim, bct
REAL(KIND(0D0)), INTENT(INOUT) :: zN
REAL(dp), INTENT(IN) :: felim, bct
REAL(dp), INTENT(INOUT) :: zN
END SUBROUTINE n2f
END INTERFACE
......@@ -147,10 +147,10 @@ CONTAINS
morfun = phy%morfun
READ (lun, NML=cmo, END=10)
10 CONTINUE
phy%A0 = A0*1D-3/cfds ! convert to m3/mmol/s
phy%A0 = A0*1E-3_dp/cfds ! convert to m3/mmol/s
phy%alpha = alpha/cfds ! convert to m2/W*mol/gChl/s
phy%Q0N = Q0N
phy%QsN = 0.5D0*Q0N
phy%QsN = 0.5E0_dp*Q0N
phy%Q0P = Q0P
phy%rdl = rdl
phy%RC = RC/cfds
......@@ -176,7 +176,7 @@ CONTAINS
CASE ('houlton')
phy%fT => ftp_houlton
CASE ('')
IF (phy%fF.EQ.0D0) THEN
IF (phy%fF.EQ.0E0_dp) THEN
phy%fT => ftp_orphy
ELSE
phy%fT => ftp_diaphy
......@@ -213,11 +213,11 @@ CONTAINS
ELSE
phy%sifun => si_tchss
END IF
IF (phy%fF.EQ.1D0) THEN
IF (phy%fF.EQ.1E0_dp) THEN
phy%fixn2 => facn2f
ELSE
phy%fixn2 => non2f
phy%FixN = 0D0
phy%FixN = 0E0_dp
END IF
IF (pf) WRITE (stdout, NML=cmo)
END SUBROUTINE read_params
......@@ -254,32 +254,32 @@ CONTAINS
V0 = phy%V0*ft*felim
V0P = phy%V0*ft ! no Fe limitation of P uptake
RC = phy%RC*ft
VNmax = V0*(1D0 - phy%Q0P/phy%QP) ! (8)
VNmax = V0*(1E0_dp - phy%Q0P/phy%QP) ! (8)
an = phy%A0*phy%DIN
ap = phy%A0*phy%DIP
! phy%VNh0 is actually vnh/vnm (to avoid division by 0 for fn)
phy%VNh0 = an/(SQRT(MAX(an, 1D-40)) + SQRT(VNmax))**2 ! (7)
phy%sVPh = SQRT(ap*V0P)/(SQRT(MAX(ap, 1D-40)) + SQRT(V0P)) ! (7)
phy%VNh0 = an/(SQRT(MAX(an, 1E-40_dp)) + SQRT(VNmax))**2 ! (7)
phy%sVPh = SQRT(ap*V0P)/(SQRT(MAX(ap, 1E-40_dp)) + SQRT(V0P)) ! (7)
zN = phy%zN
IF ((an*V0 + phy%sVPh).GT.0D0) THEN
IF ((an*V0 + phy%sVPh).GT.0E0_dp) THEN
phy%fN = phy%sVPh/(SQRT(V0*phy%VNh0**1.5*phy%Q0P/phy%QN) + phy%sVPh) ! (E42)
ELSE
phy%fN = 0.5D0
phy%fN = 0.5E0_dp
END IF
phy%fV = MAX(phy%QsN/phy%QN - zN*(phy%QN - phy%Q0N), 0D0) ! allocation for nutrient acquisition
phy%fV = MAX(phy%QsN/phy%QN - zN*(phy%QN - phy%Q0N), 0E0_dp) ! allocation for nutrient acquisition
phy%VN = phy%fV*phy%fN*VNmax*phy%VNh0
CALL phy%fixn2 (felim, bct, zN) ! diazotrophy (facn2f or non2f)
phy%fC = 1D0 - phy%fV - phy%QsN/phy%QN ! allocation for photosynthesis
phy%fC = 1E0_dp - phy%fV - phy%QsN/phy%QN ! allocation for photosynthesis
phy%mu0 = V0/(phy%rdl + phy%daylen) ! potential CO2 fixation rate
phy%I0 = RC*phy%zC/(phy%alpha*phy%daylen) ! theshold irradiance (A5)
CALL phy%sifun (RC, phin, zN) ! Chl dynamics: si_tchss (steady state) or si_dynchl (dynamic)
phy%A = phy%daylen*phy%mu0*phy%SI*(1D0 - phy%zC*phy%tch)& ! (A2)
phy%A = phy%daylen*phy%mu0*phy%SI*(1E0_dp - phy%zC*phy%tch)& ! (A2)
- RC*phy%zC*phy%tch
phy%VP = phy%fV*(1D0 - phy%fN)*phy%sVPh**2 ! (6)
phy%VP = phy%fV*(1E0_dp - phy%fN)*phy%sVPh**2 ! (6)
phy%resC = zN*(phy%VN + phy%FixN)
VC0 = phy%fC*phy%A - phy%resC ! (1)
phy%XC = MAX(VC0*phy%Q0P/phy%QP - phy%VP/phy%Q0P, 0D0)&
*MAX(2D0 - phy%QP/phy%Q0P, 0D0)
phy%XC = MAX(VC0*phy%Q0P/phy%QP - phy%VP/phy%Q0P, 0E0_dp)&
*MAX(2E0_dp - phy%QP/phy%Q0P, 0E0_dp)
! extracellular release (no DOC in UVic: simply subtract XC from VC0)
phy%VC = VC0 - phy%XC ! (A2)
END SUBROUTINE grow
......@@ -291,16 +291,16 @@ CONTAINS
REAL(dp), INTENT(INOUT) :: zN
REAL(dp) :: FNmax, fnF, fvF
FNmax = phy%F0N*felim*phy%fTNF(bct)
fnF = phy%sVPh/(SQRT(MAX(FNmax, 1D-30)*phy%Q0P/phy%QN) + phy%sVPh) ! (E54)
fvF = MAX(phy%QsN/phy%QN - phy%zF*(phy%QN - phy%Q0N), 0D0)
fnF = phy%sVPh/(SQRT(MAX(FNmax, 1E-30_dp)*phy%Q0P/phy%QN) + phy%sVPh) ! (E54)
fvF = MAX(phy%QsN/phy%QN - phy%zF*(phy%QN - phy%Q0N), 0E0_dp)
phy%FixN = fvF*fnF*(1 - phy%Q0P/phy%QP)*FNmax
IF (phy%VN.LT.phy%FixN) THEN ! fix N2
zN = phy%zF
phy%fN = fnF
phy%fV = fvF
phy%VN = 0D0
phy%VN = 0E0_dp
ELSE
phy%FixN = 0D0
phy%FixN = 0E0_dp
END IF
END SUBROUTINE facn2f
......@@ -318,15 +318,15 @@ CONTAINS
REAL(dp), INTENT(INOUT) :: RC
REAL(dp), INTENT(IN) :: phin, zN
REAL(dp) :: paravg
phy%tch = phy%theta/MAX(phy%fC, 1D-30) ! (A1)
IF (phy%PARb.GT.0D0) THEN
phy%SI = physi(phy, phin, phy%PARs, phy%PARb, 1D0)
phy%tch = phy%theta/MAX(phy%fC, 1E-30_dp) ! (A1)
IF (phy%PARb.GT.0E0_dp) THEN
phy%SI = physi(phy, phin, phy%PARs, phy%PARb, 1E0_dp)
ELSE
phy%SI = 0D0
phy%SI = 0E0_dp
END IF
paravg = (phy%PARs - phy%PARb)/phin
phy%rctht = phy%daylen*(phy%alpha*paravg*(1D0 - phy%SI)&
*(1D0/phy%zC - phy%tch) - phy%SI*phy%mu0) - RC
phy%rctht = phy%daylen*(phy%alpha*paravg*(1E0_dp - phy%SI)&
*(1E0_dp/phy%zC - phy%tch) - phy%SI*phy%mu0) - RC
phy%fcdtdq = phy%Q0N/phy%QN**2 + zN
END SUBROUTINE si_dynchl
......@@ -337,21 +337,21 @@ CONTAINS
REAL(dp), INTENT(INOUT) :: RC
REAL(dp), INTENT(IN) :: phin, zN
REAL(dp) :: paravg, PARb, aim, depfrac
IF ((phy%PARs.GT.phy%I0).AND.(phy%mu0.GT.1D-30)) THEN ! (A4)
IF ((phy%PARs.GT.phy%I0).AND.(phy%mu0.GT.1E-30_dp)) THEN ! (A4)
PARb = MAX(phy%I0, phy%PARb)
depfrac = LOG(phy%PARs/PARb)
paravg = (phy%PARs - PARb)/depfrac
aim = phy%alpha*paravg/phy%mu0
phy%tch = 1D0/phy%zC + (1D0 - lambertw((1D0 + RC/(phy%daylen*phy%mu0))&
*EXP(1D0 + MIN(aim/phy%zC, 2D2)), 0, 0))/aim
phy%tch = 1E0_dp/phy%zC + (1E0_dp - lambertw((1E0_dp + RC/(phy%daylen*phy%mu0))&
*EXP(1E0_dp + MIN(aim/phy%zC, 2E2_dp)), 0, 0))/aim
phy%SI = physi(phy, phin, phy%PARs, PARb, depfrac)
phy%theta = phy%tch*phy%fC
RC = RC*depfrac/phin ! for calculation of phy%A in SR grow
ELSE
depfrac = 0D0
phy%theta = 0D0
phy%tch = 0D0
phy%SI = 0D0
depfrac = 0E0_dp
phy%theta = 0E0_dp
phy%tch = 0E0_dp
phy%SI = 0E0_dp
END IF
END SUBROUTINE si_tchss
......@@ -363,17 +363,17 @@ CONTAINS
CLASS(ocm), INTENT(IN) :: phy
REAL(dp), INTENT(IN) :: phin, PARms, PARmb, df
REAL(dp) :: atm, SI
atm = 2D0*phy%alpha*phy%tch/phy%mu0
atm = 2E0_dp*phy%alpha*phy%tch/phy%mu0
SI = (df - (dei(-atm*PARms) - dei(-atm*PARmb)&
+ ((1D0 - EXP(-atm*PARmb))/PARmb&
- (1D0 - EXP(-atm*PARms))/PARms)/atm))/phin
+ ((1E0_dp - EXP(-atm*PARmb))/PARmb&
- (1E0_dp - EXP(-atm*PARms))/PARms)/atm))/phin
END FUNCTION physi
FUNCTION ftp_orphy (bct) RESULT (ftp)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: bct
REAL(dp) :: ftp
ftp = bct/5D0
ftp = bct/5E0_dp
END FUNCTION ftp_orphy
FUNCTION ftp_houlton (bct) RESULT (ftp)
......@@ -387,14 +387,14 @@ CONTAINS
IMPLICIT NONE
REAL(dp), INTENT(IN) :: bct
REAL(dp) :: ftp
ftp = MAX(bct - 2.6, 0D0)/2D0
ftp = MAX(bct - 2.6, 0E0_dp)/2E0_dp
END FUNCTION ftp_diaphy
! mortality as in the original UVic
SUBROUTINE mort_UVic (phyt, phyc, bct)
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(IN) :: phyc, bct
REAL(dp), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyt%mor*phyc
phyt%loss%partN = phyt%loss%partC*phyt%QN
phyt%loss%partP = phyt%loss%partC*phyt%QP
......@@ -407,7 +407,7 @@ CONTAINS
SUBROUTINE mort_D_UVic (phyt, phyc, bct)
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(IN) :: phyc, bct
REAL(dp), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyt%mort*bct*phyc
phyt%loss%partN = phyt%loss%partC*phyt%QN
phyt%loss%partP = phyt%loss%partC*phyt%QP
......@@ -420,9 +420,9 @@ CONTAINS
SUBROUTINE mort_aggreg (phyt, phyc, bct)
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyc*(phyc*phyt%mor*MAX(2D0 - phyt%QN/phyt%Q0N, 0D0)& ! aggregation
+ phyt%mort*(1D0 - EXP(-phyt%lamor*phyc))) ! sinking
REAL(dp), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyc*(phyc*phyt%mor*MAX(2E0_dp - phyt%QN/phyt%Q0N, 0E0_dp)& ! aggregation
+ phyt%mort*(1E0_dp - EXP(-phyt%lamor*phyc))) ! sinking
phyt%loss%partN = phyt%loss%partC*phyt%QN
phyt%loss%partP = phyt%loss%partC*phyt%QP
END SUBROUTINE mort_aggreg
......@@ -431,8 +431,8 @@ CONTAINS
SUBROUTINE mort_mp (phyt, phyc, bct)
IMPLICIT NONE
CLASS(ocm), INTENT(INOUT):: phyt
REAL(KIND(0D0)), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyt%mor*phyc*(1D0 - EXP(-phyt%lamor*phyc))
REAL(dp), INTENT(IN) :: phyc, bct
phyt%loss%partC = phyt%mor*phyc*(1E0_dp - EXP(-phyt%lamor*phyc))
phyt%loss%partN = phyt%loss%partC*phyt%QN
phyt%loss%partP = phyt%loss%partC*phyt%QP
END SUBROUTINE mort_mp
......
f90/lambert.f90
View file @ 749d30fc
......@@ -5,7 +5,7 @@ MODULE lambert
IMPLICIT NONE
PRIVATE
PUBLIC :: lambertw, lwm1
INTEGER, PARAMETER :: dp=