/[MITgcm]/MITgcm_contrib/bling/pkg/bling_production.F
ViewVC logotype

Contents of /MITgcm_contrib/bling/pkg/bling_production.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.5 - (show annotations) (download)
Thu May 19 18:35:05 2016 UTC (9 years, 2 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +31 -31 lines
fix calls to DIAGNOSTICS_FILL for diagnostics: BLGFERD, BLGNRD & BLGPRD
 (were 7.c instead of 8)

1 C $Header: /u/gcmpack/MITgcm_contrib/bling/pkg/bling_production.F,v 1.4 2016/05/15 00:30:35 mmazloff Exp $
2 C $Name: $
3
4 #include "BLING_OPTIONS.h"
5
6 CBOP
7 subroutine BLING_PROD(
8 I PTR_NO3, PTR_PO4, PTR_FE,
9 I PTR_O2, PTR_DON, PTR_DOP,
10 O G_NO3, G_PO4, G_FE,
11 O G_O2, G_DON, G_DOP, G_CACO3,
12 I bi, bj, imin, imax, jmin, jmax,
13 I myIter, myTime, myThid )
14
15 C =================================================================
16 C | subroutine bling_prod
17 C | o Nutrient uptake and partitioning between organic pools.
18 C | - Phytoplankton biomass-specific growth rate is calculated
19 C | as a function of light, nutrient limitation, and
20 C | temperature.
21 C | - Biomass growth xxx
22 C =================================================================
23
24 implicit none
25
26 C === Global variables ===
27 C P_sm :: Small phytoplankton biomass
28 C P_lg :: Large phytoplankton biomass
29 C P_diaz :: Diazotroph phytoplankton biomass
30
31 #include "SIZE.h"
32 #include "DYNVARS.h"
33 #include "EEPARAMS.h"
34 #include "PARAMS.h"
35 #include "GRID.h"
36 #include "BLING_VARS.h"
37 #include "PTRACERS_SIZE.h"
38 #include "PTRACERS_PARAMS.h"
39 #ifdef ALLOW_AUTODIFF
40 # include "tamc.h"
41 #endif
42
43 C === Routine arguments ===
44 C bi,bj :: tile indices
45 C iMin,iMax :: computation domain: 1rst index range
46 C jMin,jMax :: computation domain: 2nd index range
47 C myTime :: current time
48 C myIter :: current timestep
49 C myThid :: thread Id. number
50 INTEGER bi, bj, imin, imax, jmin, jmax
51 _RL myTime
52 INTEGER myIter
53 INTEGER myThid
54 C === Input ===
55 C PTR_NO3 :: nitrate concentration
56 C PTR_PO4 :: phosphate concentration
57 C PTR_FE :: iron concentration
58 C PTR_DON :: dissolved organic nitrogen concentration
59 C PTR_DOP :: dissolved organic phosphorus concentration
60 C PTR_O2 :: oxygen concentration
61 _RL PTR_NO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62 _RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63 _RL PTR_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
64 _RL PTR_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
65 _RL PTR_DON(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
66 _RL PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
67 C === Output ===
68 C G_xxx :: Tendency of xxx
69 _RL G_NO3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
70 _RL G_PO4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
71 _RL G_FE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72 _RL G_O2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73 _RL G_DON (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
74 _RL G_DOP (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75 _RL G_CACO3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
76
77 #ifdef ALLOW_BLING
78 C === Local variables ===
79 C i,j,k :: loop indicesi
80 C irr_eff :: effective irradiance
81 C NO3_lim :: nitrate limitation
82 C PO4_lim :: phosphate limitation
83 C Fe_lim :: iron limitation for phytoplankton
84 C Fe_lim_diaz :: iron limitation for diazotrophs
85 C alpha_Fe :: initial slope of the P-I curve
86 C theta_Fe :: Chl:C ratio
87 C theta_Fe_max :: Fe-replete maximum Chl:C ratio
88 C irrk :: nut-limited efficiency of algal photosystems
89 C irr_inst :: instantaneous light
90 C irr_eff :: available light
91 C mld :: mixed layer depth
92 C Pc_m :: light-saturated max photosynthesis rate for phyt
93 C Pc_m_diaz :: light-saturated max photosynthesis rate for diaz
94 C Pc_tot :: carbon-specific photosynthesis rate
95 C expkT :: temperature function
96 C mu :: net carbon-specific growth rate for phyt
97 C mu_diaz :: net carbon-specific growth rate for diaz
98 C N_uptake :: NO3 utilization by phytoplankton
99 C N_fix :: Nitrogen fixation by diazotrophs
100 C P_uptake :: PO4 utilization by phytoplankton
101 C Fe_uptake :: dissolved Fe utilization by phytoplankton
102 C CaCO3_uptake :: Calcium carbonate uptake for shell formation
103 C DON_prod :: production of dissolved organic nitrogen
104 C DOP_prod :: production of dissolved organic phosphorus
105 C O2_prod :: production of oxygen
106 C
107 INTEGER i,j,k
108 INTEGER tmp
109 _RL th1
110 _RL th2
111 _RL th3
112 _RL NO3_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113 _RL PO4_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114 _RL Fe_lim(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115 _RL Fe_lim_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116 _RL expkT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
117 _RL Pc_m
118 _RL Pc_m_diaz
119 _RL theta_Fe_max
120 _RL theta_Fe
121 _RL irrk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
122 _RL irr_inst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123 _RL irr_eff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
124 _RL mld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 _RL mu(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126 _RL mu_diaz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127 _RL PtoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128 _RL FetoN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
129 _RL N_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
130 _RL N_fix(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
131 _RL N_den_pelag(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132 _RL N_den_benthic(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
133 _RL P_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
134 _RL Fe_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
135 _RL CaCO3_uptake(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
136 _RL CaCO3_diss(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
137 _RL DON_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
138 _RL DOP_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
139 _RL DON_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
140 _RL DOP_remin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
141 _RL O2_prod(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
142 _RL frac_exp
143 _RL N_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
144 _RL P_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
145 _RL Fe_spm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
146 _RL N_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
147 _RL P_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
148 _RL Fe_dvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
149 _RL N_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
150 _RL P_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
151 _RL Fe_recycle(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
152 _RL N_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153 _RL P_reminp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
154 _RL Fe_reminsum(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
155 _RL N_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
156 _RL P_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
157 _RL Fe_remindvm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
158 _RL POC_flux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159 _RL NPP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
160 _RL NCP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
161 #ifdef ML_MEAN_PHYTO
162 _RL tmp_p_sm_ML
163 _RL tmp_p_lg_ML
164 _RL tmp_p_diaz_ML
165 _RL tmp_ML
166 #endif
167 CEOP
168
169 c ---------------------------------------------------------------------
170 c Initialize output and diagnostics
171 DO j=jmin,jmax
172 DO i=imin,imax
173 mld(i,j) = 0. _d 0
174 ENDDO
175 ENDDO
176 DO k=1,Nr
177 DO j=jmin,jmax
178 DO i=imin,imax
179 G_NO3(i,j,k) = 0. _d 0
180 G_PO4(i,j,k) = 0. _d 0
181 G_Fe(i,j,k) = 0. _d 0
182 G_O2(i,j,k) = 0. _d 0
183 G_DON(i,j,k) = 0. _d 0
184 G_DOP(i,j,k) = 0. _d 0
185 G_CaCO3(i,j,k) = 0. _d 0
186 N_uptake(i,j,k) = 0. _d 0
187 N_fix(i,j,k) = 0. _d 0
188 N_den_pelag(i,j,k) = 0. _d 0
189 N_den_benthic(i,j,k)= 0. _d 0
190 P_uptake(i,j,k) = 0. _d 0
191 Fe_uptake(i,j,k) = 0. _d 0
192 CaCO3_uptake(i,j,k) = 0. _d 0
193 DON_prod(i,j,k) = 0. _d 0
194 DOP_prod(i,j,k) = 0. _d 0
195 O2_prod(i,j,k) = 0. _d 0
196 mu_diaz(i,j,k) = 0. _d 0
197 irr_eff(i,j,k) = 0. _d 0
198 irr_inst(i,j,k) = 0. _d 0
199 irrk(i,j,k) = 0. _d 0
200 NO3_lim(i,j,k) = 0. _d 0
201 PO4_lim(i,j,k) = 0. _d 0
202 Fe_lim(i,j,k) = 0. _d 0
203 Fe_lim_diaz(i,j,k) = 0. _d 0
204 PtoN(i,j,k) = 0. _d 0
205 FetoN(i,j,k) = 0. _d 0
206 NPP(i,j,k) = 0. _d 0
207 N_reminp(i,j,k) = 0. _d 0
208 P_reminp(i,j,k) = 0. _d 0
209 Fe_reminsum(i,j,k) = 0. _d 0
210 N_remindvm(i,j,k) = 0. _d 0
211 P_remindvm(i,j,k) = 0. _d 0
212 ENDDO
213 ENDDO
214 ENDDO
215
216
217 c-----------------------------------------------------------
218 c avoid negative nutrient concentrations that can result from
219 c advection when low concentrations
220
221 #ifdef BLING_NO_NEG
222 CALL TRACER_MIN_VAL( PTR_NO3, 1. _d -7)
223 CALL TRACER_MIN_VAL( PTR_PO4, 1. _d -8)
224 CALL TRACER_MIN_VAL( PTR_FE, 1. _d -11)
225 CALL TRACER_MIN_VAL( PTR_O2, 1. _d -11)
226 CALL TRACER_MIN_VAL( PTR_DON, 1. _d -11)
227 CALL TRACER_MIN_VAL( PTR_DOP, 1. _d -11)
228 #endif
229
230 c-----------------------------------------------------------
231 c mixed layer depth calculation for light and dvm
232 c
233 CALL BLING_MIXEDLAYER(
234 U mld,
235 I bi, bj, imin, imax, jmin, jmax,
236 I myIter, myTime, myThid)
237
238 c Phytoplankton mixing
239 c The mixed layer is assumed to homogenize vertical gradients of phytoplankton.
240 c This allows for basic Sverdrup dynamics in a qualitative sense.
241 c This has not been thoroughly tested, and care should be
242 c taken with its interpretation.
243
244
245 #ifdef ML_MEAN_PHYTO
246 DO j=jmin,jmax
247 DO i=imin,imax
248
249 tmp_p_sm_ML = 0. _d 0
250 tmp_p_lg_ML = 0. _d 0
251 tmp_p_diaz_ML = 0. _d 0
252 tmp_ML = 0. _d 0
253
254 DO k=1,Nr
255
256 IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
257 IF ((-rf(k+1) .le. mld(i,j)).and.
258 & (-rf(k+1).lt.200. _d 0)) THEN
259 tmp_p_sm_ML = tmp_p_sm_ML+P_sm(i,j,k,bi,bj)*drF(k)
260 & *hFacC(i,j,k,bi,bj)
261 tmp_p_lg_ML = tmp_p_lg_ML+P_lg(i,j,k,bi,bj)*drF(k)
262 & *hFacC(i,j,k,bi,bj)
263 tmp_p_diaz_ML = tmp_p_diaz_ML+P_diaz(i,j,k,bi,bj)*drF(k)
264 & *hFacC(i,j,k,bi,bj)
265 tmp_ML = tmp_ML + drF(k)
266 ENDIF
267 ENDIF
268
269 ENDDO
270
271 DO k=1,Nr
272
273 IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
274 IF ((-rf(k+1) .le. mld(i,j)).and.
275 & (-rf(k+1).lt.200. _d 0)) THEN
276
277 P_sm(i,j,k,bi,bj) = max(1. _d -8,tmp_p_sm_ML/tmp_ML)
278 P_lg(i,j,k,bi,bj) = max(1. _d -8,tmp_p_lg_ML/tmp_ML)
279 P_diaz(i,j,k,bi,bj) = max(1. _d -8,tmp_p_diaz_ML/tmp_ML)
280
281 ENDIF
282 ENDIF
283
284 ENDDO
285 ENDDO
286 ENDDO
287
288 #endif
289
290
291 c-----------------------------------------------------------
292 c light availability for biological production
293 CALL BLING_LIGHT(
294 I mld,
295 U irr_inst, irr_eff,
296 I bi, bj, imin, imax, jmin, jmax,
297 I myIter, myTime, myThid )
298
299 c phytoplankton photoadaptation to local light level
300 DO k=1,Nr
301 DO j=jmin,jmax
302 DO i=imin,imax
303
304 irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
305 & (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
306 & min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
307
308 ENDDO
309 ENDDO
310 ENDDO
311
312 c ---------------------------------------------------------------------
313 c Nutrient uptake and partitioning between organic pools
314
315 DO k=1,Nr
316 DO j=jmin,jmax
317 DO i=imin,imax
318
319 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
320
321 c ---------------------------------------------------------------------
322 c First, calculate the limitation terms for NUT and Fe, and the
323 c Fe-limited Chl:C maximum. The light-saturated maximal photosynthesis
324 c rate term (Pc_m) is simply the product of a prescribed maximal
325 c photosynthesis rate (Pc_0), the Eppley temperature dependence, and a
326 c resource limitation term. The iron limitation term has a lower limit
327 c of Fe_lim_min and is scaled by (k_Fe2P + Fe2P_max) / Fe2P_max so that
328 c it approaches 1 as Fe approaches infinity. Thus, it is of comparable
329 c magnitude to the macro-nutrient limitation term.
330
331 c Macro-nutrient limitation
332 NO3_lim(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3)
333
334 PO4_lim(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4)
335
336 c Iron limitation
337
338 Fe_lim(i,j,k) = PTR_FE(i,j,k) / (PTR_FE(i,j,k)+k_Fe)
339
340 Fe_lim_diaz(i,j,k) = PTR_FE(i,j,k) / (PTR_FE(i,j,k)+k_Fe_diaz)
341
342 c ---------------------------------------------------------------------
343 c Diazotrophs are assumed to be more strongly temperature sensitive,
344 c given their observed restriction to relatively warm waters. Presumably
345 c this is because of the difficulty of achieving N2 fixation in an oxic
346 c environment. Thus, they have lower pc_0 and higher kappa_eppley.
347 c Taking the square root, to provide the geometric mean.
348
349 expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
350
351 c Light-saturated maximal photosynthesis rate
352
353 #ifdef BLING_ADJOINT_SAFE_tmp_xxxxxxxxxxxxxxxxxx_needs_testing
354 th1 = tanh( (NO3_lim(i,j,k)-PO4_lim(i,j,k))*1. _d 6 )
355 nut_lim = ( 1. _d 0 - th1 ) * NO3_lim(i,j,k) * 0.5 _d 0
356 & + ( 1. _d 0 + th1 ) * PO4_lim(i,j,k) * 0.5 _d 0
357
358 th2 = tanh( (nut_lim-Fe_lim(i,j,k))*1. _d 6 )
359 tot_lim = ( 1. _d 0 - th2 ) * nut_lim * 0.5 _d 0
360 & + ( 1. _d 0 + th2 ) * Fe_lim(i,j,k) * 0.5 _d 0
361
362 th3 = tanh( (PO4_lim(i,j,k)-Fe_lim(i,j,k))*1. _d 6 )
363 diaz_lim = ( 1. _d 0 - th3 ) * PO4_lim(i,j,k) * 0.5 _d 0
364 & + ( 1. _d 0 + th3 ) * Fe_lim(i,j,k) * 0.5 _d 0
365
366
367 Pc_m = Pc_0 * expkT(i,j,k) * tot_lim
368 & * maskC(i,j,k,bi,bj)
369
370 Pc_m_diaz = Pc_0_diaz
371 & * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
372 & * diaz_lim * maskC(i,j,k,bi,bj)
373
374 #else
375
376 Pc_m = Pc_0 * expkT(i,j,k)
377 & * min(NO3_lim(i,j,k), PO4_lim(i,j,k), Fe_lim(i,j,k))
378 & * maskC(i,j,k,bi,bj)
379
380 Pc_m_diaz = Pc_0_diaz
381 & * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
382 & * min(PO4_lim(i,j,k), Fe_lim_diaz(i,j,k))
383 & * maskC(i,j,k,bi,bj)
384
385 CMM( Pc_m and Pc_m_diaz crash adjoint if get too small
386 #ifdef BLING_ADJOINT_SAFE
387 Pc_m = MAX(Pc_m ,maskC(i,j,k,bi,bj)*1. _d -15)
388 Pc_m_diaz = MAX(Pc_m_diaz,maskC(i,j,k,bi,bj)*1. _d -15)
389 #endif
390 CMM)
391 #endif
392
393
394 c ---------------------------------------------------------------------
395 c Fe limitation 1) reduces photosynthetic efficiency (alpha_Fe)
396 c and 2) reduces the maximum achievable Chl:C ratio (theta_Fe)
397 c below a prescribed, Fe-replete maximum value (theta_Fe_max),
398 c to approach a prescribed minimum Chl:C (theta_Fe_min) under extreme
399 c Fe-limitation.
400
401 theta_Fe_max = theta_Fe_max_lo+
402 & (theta_Fe_max_hi-theta_Fe_max_lo)*Fe_lim(i,j,k)
403
404 theta_Fe = theta_Fe_max/(1. _d 0 + alpha_photo*theta_Fe_max
405 & *irr_mem(i,j,k,bi,bj)/(epsln + 2. _d 0*Pc_m))
406
407 c ---------------------------------------------------------------------
408 c Nutrient-limited efficiency of algal photosystems, irrk, is calculated
409 c with the iron limitation term included as a multiplier of the
410 c theta_Fe_max to represent the importance of Fe in forming chlorophyll
411 c accessory antennae, which do not affect the Chl:C but still affect the
412 c phytoplankton ability to use light (eg Stzrepek & Harrison, Nature 2004).
413
414 irrk(i,j,k) = Pc_m/(epsln + alpha_photo*theta_Fe_max) +
415 & irr_mem(i,j,k,bi,bj)/2. _d 0
416
417 c Carbon-specific photosynthesis rate
418 mu(i,j,k) = Pc_m * ( 1. _d 0 - exp(-irr_eff(i,j,k)
419 & /(epsln + irrk(i,j,k))))
420
421 mu_diaz(i,j,k) = Pc_m_diaz * ( 1. _d 0 - exp(-irr_eff(i,j,k)
422 & /(epsln + irrk(i,j,k))))
423
424 ENDIF
425 ENDDO
426 ENDDO
427 ENDDO
428
429 c Instantaneous nutrient concentration in phyto biomass
430 c Separate loop so adjoint stuff above can be outside loop
431 c (fix for recomputations)
432 CMM(
433 CADJ STORE P_sm = comlev1, key = ikey_dynamics, kind=isbyte
434 CADJ STORE P_lg = comlev1, key = ikey_dynamics, kind=isbyte
435 CADJ STORE P_diaz = comlev1, key = ikey_dynamics, kind=isbyte
436 CMM)
437 DO k=1,Nr
438 DO j=jmin,jmax
439 DO i=imin,imax
440
441 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
442
443 c expkT = exp(kappa_eppley * theta(i,j,k,bi,bj))
444
445 P_lg(i,j,k,bi,bj) = P_lg(i,j,k,bi,bj) +
446 & P_lg(i,j,k,bi,bj)*(mu(i,j,k) - lambda_0
447 & *expkT(i,j,k) *
448 & (P_lg(i,j,k,bi,bj)/pivotal)**(1. / 3.))
449 & * PTRACERS_dTLev(k)
450
451 P_sm(i,j,k,bi,bj) = P_sm(i,j,k,bi,bj) +
452 & P_sm(i,j,k,bi,bj)*(mu(i,j,k) - lambda_0
453 & *expkT(i,j,k) * (P_sm(i,j,k,bi,bj)/pivotal) )
454 & * PTRACERS_dTLev(k)
455
456 P_diaz(i,j,k,bi,bj) = P_diaz(i,j,k,bi,bj) +
457 & P_diaz(i,j,k,bi,bj)*(mu_diaz(i,j,k) - lambda_0
458 & *expkT(i,j,k) * (P_diaz(i,j,k,bi,bj)/pivotal) )
459 & * PTRACERS_dTLev(k)
460
461 ENDIF
462 ENDDO
463 ENDDO
464 ENDDO
465
466 DO k=1,Nr
467 DO j=jmin,jmax
468 DO i=imin,imax
469
470 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
471
472 c use the diagnostic biomass to calculate the chl concentration
473 chl(i,j,k,bi,bj) = max(chl_min, CtoN * 12.01 * theta_Fe *
474 & (P_lg(i,j,k,bi,bj) + P_sm(i,j,k,bi,bj)
475 & + P_diaz(i,j,k,bi,bj)))
476
477 c stoichiometry
478 PtoN(i,j,k) = PtoN_min + (PtoN_max - PtoN_min) *
479 & PTR_PO4(i,j,k) / (k_PtoN + PTR_PO4(i,j,k))
480
481 FetoN(i,j,k) = FetoN_min + (FetoN_max - FetoN_min) *
482 & PTR_FE(i,j,k) / (k_FetoN + PTR_FE(i,j,k))
483
484 c Nutrient uptake
485 N_uptake(i,j,k) = mu(i,j,k)*(P_sm(i,j,k,bi,bj)
486 & + P_lg(i,j,k,bi,bj))
487
488 N_fix(i,j,k) = mu_diaz(i,j,k) * P_diaz(i,j,k,bi,bj)
489
490 P_uptake(i,j,k) = (N_uptake(i,j,k) +
491 & N_fix(i,j,k)) * PtoN(i,j,k)
492
493 Fe_uptake(i,j,k) = (N_uptake(i,j,k) +
494 & N_fix(i,j,k)) * FetoN(i,j,k)
495
496 c ---------------------------------------------------------------------
497 c Alkalinity is consumed through the production of CaCO3. Here, this is
498 c simply a linear function of the implied growth rate of small
499 c phytoplankton, which gave a reasonably good fit to the global
500 c observational synthesis of Dunne (2009). This is consistent
501 c with the findings of Jin et al. (GBC,2006).
502
503 CaCO3_uptake(i,j,k) = P_sm(i,j,k,bi,bj) * phi_sm *expkT(i,j,k)
504 & * mu(i,j,k) * CatoN
505
506 c ---------------------------------------------------------------------
507 c Partitioning between organic pools
508
509 c The uptake of nutrients is assumed to contribute to the growth of
510 c phytoplankton, which subsequently die and are consumed by heterotrophs.
511 c This can involve the transfer of nutrient elements between many
512 c organic pools, both particulate and dissolved, with complex histories.
513 c We take a simple approach here, partitioning the total uptake into two
514 c fractions - sinking and non-sinking - as a function of temperature,
515 c following Dunne et al. (2005).
516 c Then, the non-sinking fraction is further subdivided, such that the
517 c majority is recycled instantaneously to the inorganic nutrient pool,
518 c representing the fast turnover of labile dissolved organic matter via
519 c the microbial loop, and the remainder is converted to semi-labile
520 c dissolved organic matter. Iron and macro-nutrient are treated
521 c identically for the first step, but all iron is recycled
522 c instantaneously in the second step (i.e. there is no dissolved organic
523 c iron pool).
524
525 c sinking fraction: particulate organic matter
526
527 c expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
528
529 frac_exp = (phi_sm + phi_lg * (mu(i,j,k)/
530 & (epsln + lambda_0*expkT(i,j,k)))**2.)/
531 & (1. + (mu(i,j,k)/(epsln + lambda_0*expkT(i,j,k)))**2.)*
532 & exp(kappa_remin * theta(i,j,k,bi,bj))
533
534 N_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
535 & (N_uptake(i,j,k) + N_fix(i,j,k))
536
537 P_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
538 & P_uptake(i,j,k)
539
540 Fe_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
541 & Fe_uptake(i,j,k)
542
543 N_dvm(i,j,k) = frac_exp *
544 & (N_uptake(i,j,k) + N_fix(i,j,k)) - N_spm(i,j,k)
545
546 P_dvm(i,j,k) = frac_exp * P_uptake(i,j,k) -
547 & P_spm(i,j,k)
548
549 Fe_dvm(i,j,k) = frac_exp * Fe_uptake(i,j,k) -
550 & Fe_spm(i,j,k)
551
552 c the remainder is divided between instantaneously recycled and
553 c long-lived dissolved organic matter.
554
555 DON_prod(i,j,k) = phi_DOM*(N_uptake(i,j,k)
556 & + N_fix(i,j,k) - N_spm(i,j,k)
557 & - N_dvm(i,j,k))
558
559 DOP_prod(i,j,k) = phi_DOM*(P_uptake(i,j,k)
560 & - P_spm(i,j,k) - P_dvm(i,j,k))
561
562 N_recycle(i,j,k) = N_uptake(i,j,k) + N_fix(i,j,k)
563 & - N_spm(i,j,k) - DON_prod(i,j,k)
564 & - N_dvm(i,j,k)
565
566 P_recycle(i,j,k) = P_uptake(i,j,k)
567 & - P_spm(i,j,k) - DOP_prod(i,j,k)
568 & - P_dvm(i,j,k)
569
570 Fe_recycle(i,j,k) = Fe_uptake(i,j,k)
571 & - Fe_spm(i,j,k) - Fe_dvm(i,j,k)
572
573 ENDIF
574
575 ENDDO
576 ENDDO
577 ENDDO
578
579
580 c-----------------------------------------------------------
581 c remineralization of sinking organic matter
582 CALL BLING_REMIN(
583 I PTR_NO3, PTR_FE, PTR_O2, irr_inst,
584 I N_spm, P_spm, Fe_spm, CaCO3_uptake,
585 U N_reminp, P_reminp, Fe_reminsum,
586 U N_den_benthic, CACO3_diss,
587 I bi, bj, imin, imax, jmin, jmax,
588 I myIter, myTime, myThid)
589
590 c-----------------------------------------------------------
591 c remineralization from diel vertical migration
592 CALL BLING_DVM(
593 I N_dvm,P_dvm,Fe_dvm,
594 I PTR_O2, mld,
595 O N_remindvm, P_remindvm, Fe_remindvm,
596 I bi, bj, imin, imax, jmin, jmax,
597 I myIter, myTime, myThid)
598
599
600 c-----------------------------------------------------------
601 c sub grid scale sediments
602 #ifdef USE_SGS_SED
603 CALL BLING_SGS(
604 I xxx,
605 O xxx,
606 I bi, bj, imin, imax, jmin, jmax,
607 I myIter, myTime, myThid)#endif
608 #endif
609
610
611 c-----------------------------------------------------------
612 c
613
614 DO k=1,Nr
615 DO j=jmin,jmax
616 DO i=imin,imax
617
618 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
619
620
621 c Dissolved organic matter slow remineralization
622
623 #ifdef BLING_NO_NEG
624 DON_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DON
625 & *PTR_DON(i,j,k),0. _d 0)
626 DOP_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DOP
627 & *PTR_DOP(i,j,k),0. _d 0)
628 #else
629 DON_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DON
630 & *PTR_DON(i,j,k)
631 DOP_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DOP
632 & *PTR_DOP(i,j,k)
633 #endif
634
635
636 c Pelagic denitrification
637 c If anoxic
638 cxx IF (PTR_O2(i,j,k) .lt. 0. _d 0) THEN
639
640 IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
641 IF (PTR_NO3(i,j,k) .gt. oxic_min) THEN
642 N_den_pelag(i,j,k) = max(epsln, (NO3toN *
643 & ((1. _d 0 - phi_DOM) * (N_reminp(i,j,k)
644 & + N_remindvm(i,j,k)) + DON_remin(i,j,k) +
645 & N_recycle(i,j,k))) - N_den_benthic(i,j,k))
646 ENDIF
647 ENDIF
648
649 c Carbon flux diagnostic
650 POC_flux(i,j,k) = CtoN * N_spm(i,j,k)
651
652 NPP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)) * CtoN
653
654 c oxygen production through photosynthesis
655 O2_prod(i,j,k) = O2toN * N_uptake(i,j,k)
656 & + (O2toN - 1.25 _d 0) * N_fix(i,j,k)
657
658
659
660 c-----------------------------------------------------------
661 C ADD TERMS
662
663 c Nutrients
664 c Sum of fast recycling, decay of sinking POM, and decay of DOM,
665 c less uptake, (less denitrification).
666
667 G_PO4(i,j,k) = -P_uptake(i,j,k) + P_recycle(i,j,k)
668 & + (1. _d 0 - phi_DOM) * (P_reminp(i,j,k)
669 & + P_remindvm(i,j,k)) + DOP_remin(i,j,k)
670
671 G_NO3(i,j,k) = -N_uptake(i,j,k)
672 IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
673 c Anoxic
674 G_NO3(i,j,k) = G_NO3(i,j,k)
675 & - N_den_pelag(i,j,k) - N_den_benthic(i,j,k)
676 ELSE
677 c Oxic
678 G_NO3(i,j,k) = G_NO3(i,j,k)
679 & + N_recycle(i,j,k) + (1. _d 0 - phi_DOM) *
680 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
681 & + DON_remin(i,j,k)
682 ENDIF
683
684 cxxxx check
685 NCP(i,j,k) = (-G_NO3(i,j,k) + N_fix(i,j,k)) * CtoN
686
687 c Iron
688 c remineralization, sediments and adsorption are all bundled into
689 c Fe_reminsum
690
691 G_FE(i,j,k) = -Fe_uptake(i,j,k) + Fe_reminsum(i,j,k)
692 & + Fe_remindvm(i,j,k) + Fe_recycle(i,j,k)
693
694 c Dissolved Organic Matter
695 c A fraction of POM remineralization goes into dissolved pools.
696
697 G_DON(i,j,k) = DON_prod(i,j,k) + phi_DOM *
698 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
699 & - DON_remin(i,j,k)
700
701 G_DOP(i,j,k) = DOP_prod(i,j,k) + phi_DOM *
702 & (P_reminp(i,j,k) + P_remindvm(i,j,k))
703 & - DOP_remin(i,j,k)
704
705 c Oxygen:
706 c Assuming constant O2:N ratio in terms of oxidant required per mol of organic N.
707 c This implies a constant stoichiometry of C:N and H:N (where H is reduced, organic H).
708 c Because the N provided by N2 fixation is reduced from N2, rather than NO3-, the
709 c o2_2_n_fix is slightly less than the NO3- based ratio (by 1.25 mol O2/ mol N).
710 c Account for the organic matter respired through benthic denitrification by
711 c subtracting 5/4 times the benthic denitrification NO3 utilization rate from
712 c the overall oxygen consumption.
713
714 G_O2(i,j,k) = O2_prod(i,j,k)
715 c If oxic
716 IF (PTR_O2(i,j,k) .gt. oxic_min) THEN
717 G_O2(i,j,k) = G_O2(i,j,k)
718 & -O2toN * ((1. _d 0 - phi_DOM) *
719 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
720 & + DON_remin(i,j,k) + N_recycle(i,j,k))
721 c If anoxic but NO3 concentration is very low
722 c (generate negative O2; proxy for HS-).
723 ELSEIF (PTR_NO3(i,j,k) .lt. oxic_min) THEN
724 G_O2(i,j,k) = G_O2(i,j,k)
725 & -O2toN * ((1. _d 0 - phi_DOM) *
726 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
727 & + DON_remin(i,j,k) + N_recycle(i,j,k))
728 & + N_den_benthic(i,j,k) * 1.25 _d 0
729 ENDIF
730
731 G_CaCO3(i,j,k) = CaCO3_diss(i,j,k) - CaCO3_uptake(i,j,k)
732 cxx sediments not accounted for
733
734 ENDIF
735
736 ENDDO
737 ENDDO
738 ENDDO
739
740
741 c ---------------------------------------------------------------------
742
743 #ifdef ALLOW_DIAGNOSTICS
744 IF ( useDiagnostics ) THEN
745
746 c 3d global variables
747 CALL DIAGNOSTICS_FILL(P_sm, 'BLGPSM ',0,Nr,1,bi,bj,myThid)
748 CALL DIAGNOSTICS_FILL(P_lg, 'BLGPLG ',0,Nr,1,bi,bj,myThid)
749 CALL DIAGNOSTICS_FILL(P_diaz, 'BLGPDIA ',0,Nr,1,bi,bj,myThid)
750 CALL DIAGNOSTICS_FILL(chl, 'BLGCHL ',0,Nr,1,bi,bj,myThid)
751 CALL DIAGNOSTICS_FILL(irr_mem,'BLGIMEM ',0,Nr,1,bi,bj,myThid)
752 c 3d local variables
753 CALL DIAGNOSTICS_FILL(irrk, 'BLGIRRK ',0,Nr,2,bi,bj,myThid)
754 CALL DIAGNOSTICS_FILL(irr_eff, 'BLGIEFF ',0,Nr,2,bi,bj,myThid)
755 CALL DIAGNOSTICS_FILL(Fe_lim, 'BLGFELIM',0,Nr,2,bi,bj,myThid)
756 CALL DIAGNOSTICS_FILL(NO3_lim, 'BLGNLIM ',0,Nr,2,bi,bj,myThid)
757 CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
758 CALL DIAGNOSTICS_FILL(NPP, 'BLGNPP ',0,Nr,2,bi,bj,myThid)
759 CALL DIAGNOSTICS_FILL(NCP, 'BLGNCP ',0,Nr,2,bi,bj,myThid)
760 c CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEAI',0,Nr,2,bi,bj,
761 c & myThid)
762 c CALL DIAGNOSTICS_FILL(Fe_dvm,'BLGFEDVM',0,Nr,2,bi,bj,myThid)
763 c CALL DIAGNOSTICS_FILL(Fe_sed,'BLGFESED',0,Nr,2,bi,bj,myThid)
764 CALL DIAGNOSTICS_FILL(Fe_spm, 'BLGFESPM',0,Nr,2,bi,bj,myThid)
765 CALL DIAGNOSTICS_FILL(Fe_recycle, 'BLGFEREC',0,Nr,2,bi,bj,
766 & myThid)
767 CALL DIAGNOSTICS_FILL(Fe_remindvm, 'BLGFERD ',0,Nr,2,bi,bj,
768 & myThid)
769 c CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid)
770 CALL DIAGNOSTICS_FILL(Fe_reminsum, 'BLGFEREM',0,Nr,2,bi,bj,
771 & myThid)
772 CALL DIAGNOSTICS_FILL(Fe_uptake,'BLGFEUP ',0,Nr,2,bi,bj,myThid)
773 CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
774 & myThid)
775 CALL DIAGNOSTICS_FILL(N_den_pelag, 'BLGNDENP',0,Nr,2,bi,bj,
776 & myThid)
777 CALL DIAGNOSTICS_FILL(N_dvm, 'BLGNDVM ',0,Nr,2,bi,bj,myThid)
778 CALL DIAGNOSTICS_FILL(N_fix, 'BLGNFIX ',0,Nr,2,bi,bj,myThid)
779 CALL DIAGNOSTICS_FILL(DON_prod, 'BLGDONP ',0,Nr,2,bi,bj,myThid)
780 CALL DIAGNOSTICS_FILL(N_spm, 'BLGNSPM ',0,Nr,2,bi,bj,myThid)
781 CALL DIAGNOSTICS_FILL(N_recycle,'BLGNREC ',0,Nr,2,bi,bj,myThid)
782 CALL DIAGNOSTICS_FILL(N_remindvm,'BLGNRD ',0,Nr,2,bi,bj,myThid)
783 CALL DIAGNOSTICS_FILL(N_reminp, 'BLGNREM ',0,Nr,2,bi,bj,myThid)
784 CALL DIAGNOSTICS_FILL(N_uptake, 'BLGNUP ',0,Nr,2,bi,bj,myThid)
785 CALL DIAGNOSTICS_FILL(P_dvm, 'BLGPDVM ',0,Nr,2,bi,bj,myThid)
786 CALL DIAGNOSTICS_FILL(DOP_prod, 'BLGDOPP ',0,Nr,2,bi,bj,myThid)
787 CALL DIAGNOSTICS_FILL(P_spm, 'BLGPSPM ',0,Nr,2,bi,bj,myThid)
788 CALL DIAGNOSTICS_FILL(P_recycle,'BLGPREC ',0,Nr,2,bi,bj,myThid)
789 CALL DIAGNOSTICS_FILL(P_remindvm,'BLGPRD ',0,Nr,2,bi,bj,myThid)
790 CALL DIAGNOSTICS_FILL(P_reminp, 'BLGPREM ',0,Nr,2,bi,bj,myThid)
791 CALL DIAGNOSTICS_FILL(P_uptake, 'BLGPUP ',0,Nr,2,bi,bj,myThid)
792 c CALL DIAGNOSTICS_FILL(dvm,'BLGDVM ',0,Nr,2,bi,bj,myThid)
793 CALL DIAGNOSTICS_FILL(mu, 'BLGMU ',0,Nr,2,bi,bj,myThid)
794 CALL DIAGNOSTICS_FILL(mu_diaz, 'BLGMUDIA',0,Nr,2,bi,bj,myThid)
795 c 2d local variables
796 c CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid)
797 c CALL DIAGNOSTICS_FILL(NO3_sed,'BLGNSED ',0,1,2,bi,bj,myThid)
798 c CALL DIAGNOSTICS_FILL(PO4_sed,'BLGPSED ',0,1,2,bi,bj,myThid)
799 c CALL DIAGNOSTICS_FILL(O2_sed,'BLGO2SED',0,1,2,bi,bj,myThid)
800 c these variables are currently 1d, could be 3d for diagnostics
801 c (or diag_fill could be called inside loop - which is faster?)
802 c CALL DIAGNOSTICS_FILL(frac_exp,'BLGFEXP ',0,Nr,2,bi,bj,myThid)
803 c CALL DIAGNOSTICS_FILL(irr_mix,'BLGIRRM ',0,Nr,2,bi,bj,myThid)
804 c CALL DIAGNOSTICS_FILL(irrk,'BLGIRRK ',0,Nr,2,bi,bj,myThid)
805 c CALL DIAGNOSTICS_FILL(kFe_eq_lig,'BLGPUP ',0,Nr,2,bi,bj,myThid)
806 c CALL DIAGNOSTICS_FILL(mu,'BLGMU ',0,Nr,2,bi,bj,myThid)
807 c CALL DIAGNOSTICS_FILL(mu_diaz,'BLGMUDIA',0,Nr,2,bi,bj,myThid)
808 c CALL DIAGNOSTICS_FILL(PtoN,'BLGP2N ',0,Nr,2,bi,bj,myThid)
809 c CALL DIAGNOSTICS_FILL(FetoN,'BLGFE2N ',0,Nr,2,bi,bj,myThid)
810 c CALL DIAGNOSTICS_FILL(Pc_m,'BLGPCM ',0,Nr,2,bi,bj,myThid)
811 c CALL DIAGNOSTICS_FILL(Pc_m_diaz,'BLGPCMD',0,Nr,2,bi,bj,myThid)
812 c CALL DIAGNOSTICS_FILL(theta_Fe,'BLGTHETA',0,Nr,2,bi,bj,myThid)
813 c CALL DIAGNOSTICS_FILL(theta_Fe_max,'BLGTHETM',0,Nr,2,bi,bj,myThid)
814 c CALL DIAGNOSTICS_FILL(wsink,'BLGWSINK',0,Nr,2,bi,bj,myThid)
815 c CALL DIAGNOSTICS_FILL(zremin,'BLGZREM ',0,Nr,2,bi,bj,myThid)
816 c CALL DIAGNOSTICS_FILL(z_dvm,'BLGZDVM ',0,Nr,2,bi,bj,myThid)
817
818 ENDIF
819 #endif /* ALLOW_DIAGNOSTICS */
820
821 #endif /* ALLOW_BLING */
822
823 RETURN
824
825 END

  ViewVC Help
Powered by ViewVC 1.1.22