-
Notifications
You must be signed in to change notification settings - Fork 0
/
PhytoWater.f90
295 lines (270 loc) · 27.2 KB
/
PhytoWater.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
Subroutine PhytoWater(HydroParam,MeshParam,MeteoParam,LimnoParam,dt,dtday)
Use MeshVars
Use Hydrodynamic
Use LimnologyVars
Use Meteorological
Use uTempFunction
Implicit None
type(MeshGridParam) :: MeshParam
type(MeteorologicalParam) :: MeteoParam
type(HydrodynamicParam) :: HydroParam
type(LimnologyParam) :: LimnoParam
! Declaracao das variaveis da subrotina
Integer:: index,iElem,iLayer,gg
Real:: wDSetPhyt_Layer,wCSetPhyt_Layer,wNSetPhyt_Layer,wPSetPhyt_Layer,tDResusPhyt_bottom,tCResusPhyt_bottom,tNResusPhyt_bottom,tPResusPhyt_bottom
Real:: V
Real:: dt,dtday
Real:: NearZero = 1e-10
Do gg = 1, LimnoParam%nphy
Do iElem = 1,MeshParam%nElem
Do iLayer = HydroParam%ElCapitalM(iElem),HydroParam%ElSmallm(iElem),-1
! Light
If ( HydroParam%ElSmallm(iElem) == HydroParam%ElCapitalM(iElem) ) Then ! Only One Vertical Layer
! light_at_the_surface,in FABM it's at the top of the current layer (W/m² PAR)
LimnoParam%uLPAR0(iLayer) = MeshParam%CREDV(iElem)*(1.-HydroParam%ALB)*LimnoParam%fPAR*MeteoParam%SolarRad(iElem)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*(HydroParam%Ze(HydroParam%ElCapitalM(iElem)+1,iElem) - HydroParam%Ze(iLayer+1,iElem)))
! light_at_the_bottom,in FABM it's at the bottom of the current layer (W/m² PAR)
LimnoParam%aLPARBot(iLayer) = LimnoParam%uLPAR0(iLayer)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*1.)
Else
If ( iLayer == HydroParam%ElSmallm(iElem) ) Then !Bottom layer
! light_at_the_surface,in FABM it's at the top of the current layer (W/m² PAR)
LimnoParam%uLPAR0(iLayer) = LimnoParam%aLPARBot(iLayer+1)
! light_at_the_bottom,in FABM it's at the bottom of the current layer (W/m² PAR)
LimnoParam%aLPARBot(iLayer) = LimnoParam%uLPAR0(iLayer)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*(HydroParam%Ze(HydroParam%ElCapitalM(iElem)+1,iElem) - HydroParam%Ze(iLayer,iElem)))
Elseif ( iLayer == HydroParam%ElCapitalM(iElem) ) Then !Surface layer
! light_at_the_surface,in FABM it's at the top of the current layer (W/m² PAR)
LimnoParam%uLPAR0(iLayer) = MeshParam%CREDV(iElem)*(1.-HydroParam%ALB)*LimnoParam%fPAR*MeteoParam%SolarRad(iElem)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*(HydroParam%Ze(HydroParam%ElCapitalM(iElem)+1,iElem) - HydroParam%Ze(iLayer+1,iElem)))
! light_at_the_bottom,in FABM it's at the bottom of the current layer (W/m² PAR)
LimnoParam%aLPARBot(iLayer) = LimnoParam%uLPAR0(iLayer)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*Max(1.,HydroParam%Ze(HydroParam%ElCapitalM(iElem)+1,iElem) - HydroParam%Ze(iLayer,iElem)))
Else !Intermediary layer
! light_at_the_surface,in FABM it's at the top of the current layer (W/m² PAR)
LimnoParam%uLPAR0(iLayer) = LimnoParam%aLPARBot(iLayer+1)
! light_at_the_bottom,in FABM it's at the bottom of the current layer (W/m² PAR)
LimnoParam%aLPARBot(iLayer) = LimnoParam%uLPAR0(iLayer)*exp(-LimnoParam%aExtCoef(iLayer,iElem)*(HydroParam%Ze(HydroParam%ElCapitalM(iElem)+1,iElem) - HydroParam%Ze(iLayer,iElem)))
EndIf
Endif
EndDo
Do iLayer = HydroParam%ElSmallm(iElem), HydroParam%ElCapitalM(iElem)
!-------------------------------------------------------------------------------------!
! Phytoplankton Processes in Water !
!-------------------------------------------------------------------------------------!
!-Phyto-------------------------------------------------------------------------------!
! 1. Primary Production/Nutrient Uptake !
! 2. Respiration/Excretion !
! 3. Lise (Mortality) !
! 4. Sedimentation !
! 5. Resuspension !
!-------------------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
!1) Primary Production
! Temperature function for pelagic phytoplankton
LimnoParam%uFunTmPhyt = uFunTmBio(LimnoParam%sDTempW(iLayer,iElem),LimnoParam%cSigTmPhyt(gg),LimnoParam%cTmOptPhyt(gg))
!-----------------------------------------------------------------------
! Light functions for pelagic phytoplankton
!-----------------------------------------------------------------------
! light limitation function for each group, each group has two different types
If (LimnoParam%LightMethodPhyt == 0) Then ! case 1, without photoinhibition,Chalker (1980) model
! half-satuation light intensity at current temperature
LimnoParam%uhLPhyt = LimnoParam%hLRefPhyt(gg) * LimnoParam%uFunTmPhyt
! light limitation function for green algae, no photo-inhibition
If (HydroParam%ElSmallm(iElem)== HydroParam%ElCapitalM(iElem)) Then ! 2D - Production in the first meter
LimnoParam%aLLimPhyt = 1.0 / (LimnoParam%aExtCoef(iLayer,iElem) * 1.) * log((1.0 + LimnoParam%uLPAR0(iLayer) / (LimnoParam%uhLPhyt+NearZero)) / (1.0 + LimnoParam%aLPARBot(iLayer) / (LimnoParam%uhLPhyt+NearZero))) !função Lehman para luz (integração)
Else
LimnoParam%aLLimPhyt = 1.0 / (LimnoParam%aExtCoef(iLayer,iElem) * Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem))) * log((1.0 + LimnoParam%uLPAR0(iLayer) / (LimnoParam%uhLPhyt+NearZero)) / (1.0 + LimnoParam%aLPARBot(iLayer) / (LimnoParam%uhLPhyt+NearZero))) !função Lehman para luz (integração)
EndIf
ElseIf (LimnoParam%LightMethodPhyt == 1) Then ! case 2, Klepper et al. (1988) / Ebenhoh et al. (1997) model.
LimnoParam%uhLPhyt = LimnoParam%cLOptRefPhyt(gg) * LimnoParam%uFunTmPhyt
If (HydroParam%ElSmallm(iElem)== HydroParam%ElCapitalM(iElem)) Then ! 2D - Production in the first meter
LimnoParam%aLLimPhyt = exp(1.0) /(LimnoParam%aExtCoef(iLayer,iElem) * 1.) *(exp(- LimnoParam%aLPARBot(iLayer) /(LimnoParam%uhLPhyt+NearZero)) - exp(- LimnoParam%uLPAR0(iLayer) /(LimnoParam%uhLPhyt+NearZero)))
Else
LimnoParam%aLLimPhyt = exp(1.0) /(LimnoParam%aExtCoef(iLayer,iElem) * Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem))) *(exp(- LimnoParam%aLPARBot(iLayer) /(LimnoParam%uhLPhyt+NearZero)) - exp(- LimnoParam%uLPAR0(iLayer) /(LimnoParam%uhLPhyt+NearZero)))
EndIf
Else
Stop 'Light limitation function for Diatoms does not exist'
EndIf
! growth_rate_at_current_light_AND_temp.
LimnoParam%aMuTmLPhyt = LimnoParam%ufDay *(1.0 - LimnoParam%afCovSurfVeg)* LimnoParam%aLLimPhyt * LimnoParam%uFunTmPhyt*LimnoParam%cMuMaxPhyt(gg)
!-----------------------------------------------------------------------
! Nutrient limitation functions
!-----------------------------------------------------------------------
! Nutrient functions for diatom
! Droop_function(C)_for_Phytoplankton
LimnoParam%aCLimPhyt = max(0.0,(1.0 - LimnoParam%cCDPhytMin(gg) / (LimnoParam%rCDPhytW(iLayer,iElem,gg) + NearZero)) * LimnoParam%cCDPhytMax(gg) /(LimnoParam%cCDPhytMax(gg) - LimnoParam%cCDPhytMin(gg)))
! Droop_function(P)_for_Phytoplankton
LimnoParam%aPLimPhyt = max(0.0,(1.0 - LimnoParam%cPDPhytMin(gg) / (LimnoParam%rPDPhytW(iLayer,iElem,gg) + NearZero)) * LimnoParam%cPDPhytMax(gg) /(LimnoParam%cPDPhytMax(gg) - LimnoParam%cPDPhytMin(gg)))
! Droop_function(N)_for_Phytoplankton
LimnoParam%aNLimPhyt = max(0.0,(1.0 - LimnoParam%cNDPhytMin(gg) / (LimnoParam%rNDPhytW(iLayer,iElem,gg) + NearZero)) * LimnoParam%cNDPhytMax(gg) /(LimnoParam%cNDPhytMax(gg) - LimnoParam%cNDPhytMin(gg)))
! silica_dependence_of_growth_rate
LimnoParam%aSiLimPhyt = LimnoParam%sSiO2W(iLayer,iElem) /(LimnoParam%hSiAssDiat(gg) + LimnoParam%sSiO2W(iLayer,iElem))
! nutrient_limitation_function_of_Diatom
! aNutLimDiat = min(aPLimDiat,aNLimDiat,aSiLimDiat) ! v5.09
LimnoParam%aNutLimPhyt = min(LimnoParam%aPLimPhyt,LimnoParam%aNLimPhyt) ! pl613
!-----------------------------------------------------------------------
! Algae growth_DW
!-----------------------------------------------------------------------
! growth_rate_Phytoplankton
LimnoParam%aMuPhyt = LimnoParam%aNutLimPhyt * LimnoParam%aMuTmLPhyt
! assimilation_Algae
LimnoParam%wDAssPhyt(iLayer,iElem,gg) = LimnoParam%aMuPhyt*LimnoParam%sDPhytW(iLayer,iElem,gg)
!-------------------------------------------------------------------------------------------------------------
! Algae uptake_C
!-------------------------------------------------------------------------------------------------------------
! maximum_P_uptake_rate_of_Algae,corrected_for_P/D_ratio
LimnoParam%aVCUptMaxCorPhyt = max(0.0,LimnoParam%cVCUptMaxPhyt(gg) * LimnoParam%aLLimPhyt * LimnoParam%uFunTmPhyt *(LimnoParam%cCDPhytMax(gg) - LimnoParam%rCDPhytW(iLayer,iElem,gg))&
&/(LimnoParam%cCDPhytMax(gg) - LimnoParam%cCDPhytMin(gg)))
! C_uptake_rate_of_Algae
LimnoParam%aVCUptPhyt = LimnoParam%sDICW(iLayer,iElem) * LimnoParam%aVCUptMaxCorPhyt /(LimnoParam%cVCUptMaxPhyt(gg) / LimnoParam%cAffCUptPhyt(gg) + LimnoParam%sDICW(iLayer,iElem))
! C_uptake_Algae
LimnoParam%wCUptPhyt(iLayer,iElem,gg) = LimnoParam%aVCUptPhyt * LimnoParam%sDPhytW(iLayer,iElem,gg)
!-------------------------------------------------------------------------------------------------------------
! Algae uptake_P
!-------------------------------------------------------------------------------------------------------------
! maximum_P_uptake_rate_of_Algae,corrected_for_P/D_ratio
LimnoParam%aVPUptMaxCorPhyt = max(0.0,LimnoParam%cVPUptMaxPhyt(gg) * LimnoParam%uFunTmPhyt *(LimnoParam%cPDPhytMax(gg) - LimnoParam%rPDPhytW(iLayer,iElem,gg))&
&/(LimnoParam%cPDPhytMax(gg) - LimnoParam%cPDPhytMin(gg)))
! P_uptake_rate_of_Algae
LimnoParam%aVPUptPhyt = LimnoParam%sPO4W(iLayer,iElem) * LimnoParam%aVPUptMaxCorPhyt /(LimnoParam%cVPUptMaxPhyt(gg) / LimnoParam%cAffPUptPhyt(gg) + LimnoParam%sPO4W(iLayer,iElem))
! P_uptake_Algae
LimnoParam%wPUptPhyt(iLayer,iElem,gg) = LimnoParam%aVPUptPhyt * LimnoParam%sDPhytW(iLayer,iElem,gg)
!-------------------------------------------------------------------------------------------------------------
! Algae uptake_N
!-------------------------------------------------------------------------------------------------------------
! maximum_N_uptake_rate_of_Algae,corrected_for_N/D_ratio
LimnoParam%aVNUptMaxCorPhyt = max(0.0,LimnoParam%cVNUptMaxPhyt(gg) * LimnoParam%uFunTmPhyt *(LimnoParam%cNDPhytMax(gg) - LimnoParam%rNDPhytW(iLayer,iElem,gg))&
&/(LimnoParam%cNDPhytMax(gg) - LimnoParam%cNDPhytMin(gg)))
! half-sat._NDissW_for_uptake_by_Algae
LimnoParam%ahNUptPhyt = LimnoParam%aVNUptMaxCorPhyt / LimnoParam%cAffNUptPhyt(gg)
! N_uptake_rate_of_Algae
LimnoParam%aVNUptPhyt = (LimnoParam%sNO3W(iLayer,iElem) + LimnoParam%sNH4W(iLayer,iElem)) * LimnoParam%aVNUptMaxCorPhyt /(LimnoParam%ahNUptPhyt + LimnoParam%sNO3W(iLayer,iElem) + LimnoParam%sNH4W(iLayer,iElem))
! N_uptake_Algae
LimnoParam%wNUptPhyt(iLayer,iElem,gg) = LimnoParam%aVNUptPhyt * LimnoParam%sDPhytW(iLayer,iElem,gg)
!-----------------------------------------------------------------------
! NH4 exchange with abiotic module
!-----------------------------------------------------------------------
! fraction_ammonium_uptake_by_Algae
LimnoParam%afNH4UptPhyt = LimnoParam%sNH4W(iLayer,iElem) * LimnoParam%sNO3W(iLayer,iElem) /((LimnoParam%ahNUptPhyt + LimnoParam%sNH4W(iLayer,iElem)) *(LimnoParam%ahNUptPhyt + LimnoParam%sNO3W(iLayer,iElem))) + LimnoParam%sNH4W(iLayer,iElem) &
& * LimnoParam%ahNUptPhyt /((LimnoParam%sNH4W(iLayer,iElem) + LimnoParam%sNO3W(iLayer,iElem)) *(LimnoParam%ahNUptPhyt + LimnoParam%sNO3W(iLayer,iElem)))
! ammonium_uptake_by_Algae
LimnoParam%wNUptNH4Phyt(iLayer,iElem,gg) = LimnoParam%afNH4UptPhyt * LimnoParam%wNUptPhyt(iLayer,iElem,gg)
!-----------------------------------------------------------------------
! NO3 exchange with abiotic module
!-----------------------------------------------------------------------
! nitrate_uptake_by_Algae
LimnoParam%wNUptNO3Phyt(iLayer,iElem,gg) = LimnoParam%wNUptPhyt(iLayer,iElem,gg) - LimnoParam%wNUptNH4Phyt(iLayer,iElem,gg)
!2) Algae respiration_DW and excretion of nutrients
! temp._corrected_respiration_constant_of_Algae
LimnoParam%ukDRespbTmPhyt = LimnoParam%kDRespbPhyt(gg) * LimnoParam%uFunTmPhyt
! basal respiration_of_Algae_in_water
LimnoParam%wDRespbPhyt = LimnoParam%ukDRespbTmPhyt * LimnoParam%sDPhytW(iLayer,iElem,gg) ! Used in O2 exchange
! Photorespiration_of_Algae_in_water (ERSEM, Polimene et al.(2006)) - alfap -> fracao da producao primaria bruta (PPG)
LimnoParam%wDResppPhyt = LimnoParam%alfap(gg) * LimnoParam%wDAssPhyt(iLayer,iElem,gg)
! Total respiration_of_Algae_in_water
LimnoParam%wDRespPhyt(iLayer,iElem,gg) = LimnoParam%wDRespbPhyt + LimnoParam%wDResppPhyt
!DIC Respiration
! C_excretion_Algae_in_water (Exudation)
! wPExcrDiatW = rPDDiatW /(self%cPDDiatMin + rPDDiatW) * rPDDiatW * wDRespDiatW ! v5.09
LimnoParam%wCExcrPhytW(iLayer,iElem,gg) = (2.0 * LimnoParam%rCDPhytW(iLayer,iElem,gg)) /(LimnoParam%cCDPhytMax(gg) + LimnoParam%rCDPhytW(iLayer,iElem,gg)) * LimnoParam%rCDPhytW(iLayer,iElem,gg) * LimnoParam%wDRespPhyt(iLayer,iElem,gg)
! P_Exudation_Algae_in_water
! wPExcrDiatW = rPDDiatW /(self%cPDDiatMin + rPDDiatW) * rPDDiatW * wDRespDiatW ! v5.09
LimnoParam%wPExcrPhytW(iLayer,iElem,gg) = (2.0 * LimnoParam%rPDPhytW(iLayer,iElem,gg)) /(LimnoParam%cPDPhytMax(gg) + LimnoParam%rPDPhytW(iLayer,iElem,gg)) * LimnoParam%rPDPhytW(iLayer,iElem,gg) * LimnoParam%wDRespPhyt(iLayer,iElem,gg)
! N_Exudation_Algae_in_water
! wNExcrDiatW = rNDDiatW /(self%cNDDiatMin + rNDDiatW) * rNDDiatW * wDRespDiatW ! V5.09
LimnoParam%wNExcrPhytW(iLayer,iElem,gg) = (2.0 * LimnoParam%rNDPhytW(iLayer,iElem,gg)) /(LimnoParam%cNDPhytMax(gg) + LimnoParam%rNDPhytW(iLayer,iElem,gg)) * LimnoParam%rNDPhytW(iLayer,iElem,gg) * LimnoParam%wDRespPhyt(iLayer,iElem,gg) ! pl613
!3) Lise_Algae_in_water
LimnoParam%wDLisPhyt(iLayer,iElem,gg) = LimnoParam%cLisPhytW(gg) * LimnoParam%sDPhytW(iLayer,iElem,gg)
LimnoParam%wCLisPhyt(iLayer,iElem,gg) = LimnoParam%wDLisPhyt(iLayer,iElem,gg) * LimnoParam%rCDPhytW(iLayer,iElem,gg)
LimnoParam%wNLisPhyt(iLayer,iElem,gg) = LimnoParam%wDLisPhyt(iLayer,iElem,gg) * LimnoParam%rNDPhytW(iLayer,iElem,gg)
LimnoParam%wPLisPhyt(iLayer,iElem,gg) = LimnoParam%wDLisPhyt(iLayer,iElem,gg) * LimnoParam%rPDPhytW(iLayer,iElem,gg)
If (LimnoParam%aphy(gg)==LimnoParam%MDIAT.or.LimnoParam%aphy(gg)==LimnoParam%FDIAT) Then
LimnoParam%wSiLisPhyt(iLayer,iElem,gg) = LimnoParam%wDLisPhyt(iLayer,iElem,gg) * LimnoParam%rSiDPhytW(iLayer,iElem,gg)
Else
LimnoParam%wSiLisPhyt(iLayer,iElem,gg) = 0.
EndIf
!4) Algae Sedimentation
If ( iLayer == HydroParam%ElCapitalM(iElem) ) Then ! Bottom layer or intermediary layer has contribuition of the up layer
wDSetPhyt_Layer = LimnoParam%tDSetPhyt(iLayer,iElem,gg)/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem))
wCSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rCDPhytW(iLayer,iElem,gg)
wNSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rNDPhytW(iLayer,iElem,gg)
wPSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rPDPhytW(iLayer,iElem,gg)
Else
wDSetPhyt_Layer = LimnoParam%tDSetPhyt(iLayer,iElem,gg)/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem)) - LimnoParam%tDSetPhyt(iLayer+1,iElem,gg)/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer+1,iElem))
wCSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rCDPhytW(iLayer,iElem,gg)
wNSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rNDPhytW(iLayer,iElem,gg)
wPSetPhyt_Layer = wDSetPhyt_Layer * LimnoParam%rPDPhytW(iLayer,iElem,gg)
Endif
!5) Algae Resuspension
If (iLayer == HydroParam%ElSmallm(iElem)) Then
tDResusPhyt_bottom = LimnoParam%tDResusPhyt(iElem,gg)
tCResusPhyt_bottom = tDResusPhyt_bottom * LimnoParam%rCDPhytW(iLayer,iElem,gg)
tNResusPhyt_bottom = tDResusPhyt_bottom * LimnoParam%rNDPhytW(iLayer,iElem,gg)
tPResusPhyt_bottom = tDResusPhyt_bottom * LimnoParam%rPDPhytW(iLayer,iElem,gg)
Else
tDResusPhyt_bottom = 0.
tCResusPhyt_bottom = 0.
tNResusPhyt_bottom = 0.
tPResusPhyt_bottom = 0.
EndIf
!-------------------------------------------------------------------------------------!
! Source/Sink Term for Algae !
!-------------------------------------------------------------------------------------!
HydroParam%dVarEst(iLayer,iElem,1) = LimnoParam%wDAssPhyt(iLayer,iElem,gg) & !Assimilation
- LimnoParam%wDRespPhyt(iLayer,iElem,gg) & !Respiration
- LimnoParam%wDLisPhyt(iLayer,iElem,gg) & !Lise Celular
- wDSetPhyt_Layer & !Sedimentation
+ tDResusPhyt_bottom/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem)) & !Resuspension
- SUM(LimnoParam%wDGrzZooPhyt(iLayer,iElem,:,gg)) & !Grazing by Zooplankton
- SUM(LimnoParam%wDGrzFiAdPhyt(iLayer,iElem,:,gg)) & !Grazing by Aduld Fish
- SUM(LimnoParam%wDGrzFiAdPhyt(iLayer,iElem,:,gg)) !Grazing by Juveline Fish
HydroParam%dVarEst(iLayer,iElem,2) = LimnoParam%wCUptPhyt(iLayer,iElem,gg) & !Nutrient Uptake
- LimnoParam%wCExcrPhytW(iLayer,iElem,gg) & !Excretion (Exudation)
- LimnoParam%wCLisPhyt(iLayer,iElem,gg) & !Lise Celular
- wCSetPhyt_Layer & !Sedimentation
+ tCResusPhyt_bottom/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem)) & !Resuspension
- SUM(LimnoParam%wCGrzZooPhyt(iLayer,iElem,:,gg)) & !Grazing by Zooplankton
- SUM(LimnoParam%wCGrzFiAdPhyt(iLayer,iElem,:,gg)) & !Grazing by Aduld Fish
- SUM(LimnoParam%wCGrzFiAdPhyt(iLayer,iElem,:,gg)) !Grazing by Juveline Fish
HydroParam%dVarEst(iLayer,iElem,3) = LimnoParam%wNUptPhyt(iLayer,iElem,gg) & !Nutrient Uptake
- LimnoParam%wNExcrPhytW(iLayer,iElem,gg) & !Excretion (Exudation)
- LimnoParam%wNLisPhyt(iLayer,iElem,gg) & !Lise Celular
- wNSetPhyt_Layer & !Sedimentation
+ tNResusPhyt_bottom/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem)) & !Resuspension
- SUM(LimnoParam%wNGrzZooPhyt(iLayer,iElem,:,gg)) & !Grazing by Zooplankton
- SUM(LimnoParam%wNGrzFiAdPhyt(iLayer,iElem,:,gg)) & !Grazing by Aduld Fish
- SUM(LimnoParam%wNGrzFiAdPhyt(iLayer,iElem,:,gg)) !Grazing by Juveline Fish
HydroParam%dVarEst(iLayer,iElem,4) = LimnoParam%wPUptPhyt(iLayer,iElem,gg) & !Nutrient Uptake
- LimnoParam%wPExcrPhytW(iLayer,iElem,gg) & !Excretion (Exudation)
- LimnoParam%wPLisPhyt(iLayer,iElem,gg) & !Lise Celular
- wPSetPhyt_Layer & !Sedimentation
+ tPResusPhyt_bottom/Max(HydroParam%Pcri,HydroParam%Dzi(iLayer,iElem)) & !Resuspension
- SUM(LimnoParam%wPGrzZooPhyt(iLayer,iElem,:,gg)) & !Grazing by Zooplankton
- SUM(LimnoParam%wPGrzFiAdPhyt(iLayer,iElem,:,gg)) & !Grazing by Aduld Fish
- SUM(LimnoParam%wPGrzFiAdPhyt(iLayer,iElem,:,gg)) !Grazing by Juveline Fish
EndDo !Loop Layer
EndDo !Loop Cell
!-------------------------------------------------------------------------------------!
! Solver Transport Equation for Phytoplankton !
!-------------------------------------------------------------------------------------!
index = 0
If (LimnoParam%iTranspFlag == 0.and.LimnoParam%iAdaptativeTimeStep==0) Then ! UpWind CWC Scheme
Call UpWindCWC(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,1),LimnoParam%sDPhytW(:,:,gg),LimnoParam%sDPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,2),LimnoParam%sCPhytW(:,:,gg),LimnoParam%sCPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,3),LimnoParam%sNPhytW(:,:,gg),LimnoParam%sNPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,4),LimnoParam%sPPhytW(:,:,gg),LimnoParam%sPPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
ElseIf (LimnoParam%iTranspFlag > 0.and.LimnoParam%iAdaptativeTimeStep==0) Then ! UpWind CWC High Resolution Scheme
Call UpWindCWC_HR(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,1),LimnoParam%sDPhytW(:,:,gg),LimnoParam%sDPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,2),LimnoParam%sCPhytW(:,:,gg),LimnoParam%sCPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,3),LimnoParam%sNPhytW(:,:,gg),LimnoParam%sNPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,4),LimnoParam%sPPhytW(:,:,gg),LimnoParam%sPPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
ElseIf (LimnoParam%iTranspFlag == 0.and.LimnoParam%iAdaptativeTimeStep==1) Then ! UpWind CWC with Global Time Stepping Scheme
Call UpWindCWC_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,1),LimnoParam%sDPhytW(:,:,gg),LimnoParam%sDPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,2),LimnoParam%sCPhytW(:,:,gg),LimnoParam%sCPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,3),LimnoParam%sNPhytW(:,:,gg),LimnoParam%sNPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
Call UpWindCWC_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,4),LimnoParam%sPPhytW(:,:,gg),LimnoParam%sPPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam)
ElseIf (LimnoParam%iTranspFlag > 0.and.LimnoParam%iAdaptativeTimeStep==1) Then ! UpWind CWC High Resolution with Global Time Stepping Scheme
Call UpWindCWC_HR_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,1),LimnoParam%sDPhytW(:,:,gg),LimnoParam%sDPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,2),LimnoParam%sCPhytW(:,:,gg),LimnoParam%sCPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,3),LimnoParam%sNPhytW(:,:,gg),LimnoParam%sNPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
Call UpWindCWC_HR_GATS(index,HydroParam%uLoadVarEst,HydroParam%dVarEst(:,:,4),LimnoParam%sPPhytW(:,:,gg),LimnoParam%sPPhytWP(:,:,gg),dt,dtday,HydroParam,MeshParam,LimnoParam%iTranspFlag)
EndiF
EndDo !Loop Functional Group
Return
End