/[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.2 - (show annotations) (download)
Sun Feb 28 21:49:24 2016 UTC (9 years, 5 months ago) by mmazloff
Branch: MAIN
Changes since 1.1: +662 -144 lines
Update to BLING version 2

1 C $Header: $
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 PtoN(i,j,k) = 0. _d 0
200 FetoN(i,j,k) = 0. _d 0
201 NPP(i,j,k) = 0. _d 0
202 N_reminp(i,j,k) = 0. _d 0
203 P_reminp(i,j,k) = 0. _d 0
204 Fe_reminsum(i,j,k) = 0. _d 0
205 N_remindvm(i,j,k) = 0. _d 0
206 P_remindvm(i,j,k) = 0. _d 0
207 ENDDO
208 ENDDO
209 ENDDO
210
211
212 c-----------------------------------------------------------
213 c avoid negative nutrient concentrations that can result from
214 c advection when low concentrations
215
216 #ifdef BLING_NO_NEG
217 CALL TRACER_MIN_VAL( PTR_NO3, 1. _d -7)
218 CALL TRACER_MIN_VAL( PTR_PO4, 1. _d -8)
219 CALL TRACER_MIN_VAL( PTR_FE, 1. _d -11)
220 #endif
221
222
223 c-----------------------------------------------------------
224 c mixed layer depth calculation for light and dvm
225 c
226 CALL BLING_MIXEDLAYER(
227 U mld,
228 I bi, bj, imin, imax, jmin, jmax,
229 I myIter, myTime, myThid)
230
231
232 c Phytoplankton mixing
233 c The mixed layer is assumed to homogenize vertical gradients of phytoplankton.
234 c This allows for basic Sverdrup dynamics in a qualitative sense.
235 c This has not been thoroughly tested, and care should be
236 c taken with its interpretation.
237
238 #ifdef ML_MEAN_PHYTO
239 DO j=jmin,jmax
240 DO i=imin,imax
241
242 tmp_p_sm_ML = 0. _d 0
243 tmp_p_lg_ML = 0. _d 0
244 tmp_p_diaz_ML = 0. _d 0
245 tmp_ML = 0. _d 0
246
247 DO k=1,Nr
248
249 IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
250 IF ((-rf(k+1) .le. mld(i,j)).and.
251 & (-rf(k+1).lt.200. _d 0)) THEN
252 tmp_p_sm_ML = tmp_p_sm_ML+P_sm(i,j,k,bi,bj)*drF(k)
253 & *hFacC(i,j,k,bi,bj)
254 tmp_p_lg_ML = tmp_p_lg_ML+P_lg(i,j,k,bi,bj)*drF(k)
255 & *hFacC(i,j,k,bi,bj)
256 tmp_p_diaz_ML = tmp_p_diaz_ML+P_diaz(i,j,k,bi,bj)*drF(k)
257 & *hFacC(i,j,k,bi,bj)
258 tmp_ML = tmp_ML + drF(k)
259 ENDIF
260 ENDIF
261
262 ENDDO
263
264 DO k=1,Nr
265
266 IF (hFacC(i,j,k,bi,bj).gt.0. _d 0) THEN
267 IF ((-rf(k+1) .le. mld(i,j)).and.
268 & (-rf(k+1).lt.200. _d 0)) THEN
269
270 P_sm(i,j,k,bi,bj) = max(1. _d -8,tmp_p_sm_ML/tmp_ML)
271 P_lg(i,j,k,bi,bj) = max(1. _d -8,tmp_p_lg_ML/tmp_ML)
272 P_diaz(i,j,k,bi,bj) = max(1. _d -8,tmp_p_diaz_ML/tmp_ML)
273
274 ENDIF
275 ENDIF
276
277 ENDDO
278 ENDDO
279 ENDDO
280
281 #endif
282
283
284 c-----------------------------------------------------------
285 c light availability for biological production
286 CALL BLING_LIGHT(
287 I mld,
288 U irr_inst, irr_eff,
289 I bi, bj, imin, imax, jmin, jmax,
290 I myIter, myTime, myThid )
291
292
293
294 c phytoplankton photoadaptation to local light level
295 DO k=1,Nr
296 DO j=jmin,jmax
297 DO i=imin,imax
298
299 irr_mem(i,j,k,bi,bj) = irr_mem(i,j,k,bi,bj) +
300 & (irr_eff(i,j,k) - irr_mem(i,j,k,bi,bj))*
301 & min( 1. _d 0, gamma_irr_mem*PTRACERS_dTLev(k) )
302
303 ENDDO
304 ENDDO
305 ENDDO
306
307
308 c ---------------------------------------------------------------------
309 c Nutrient uptake and partitioning between organic pools
310
311 C!! needed??
312 C$TAF STORE P_sm = comlev1, key = ikey_dynamics, kind=isbyte
313 C$TAF STORE P_lg = comlev1, key = ikey_dynamics, kind=isbyte
314 C$TAF STORE P_diaz = comlev1, key = ikey_dynamics, kind=isbyte
315
316 DO k=1,Nr
317 DO j=jmin,jmax
318 DO i=imin,imax
319
320 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
321
322 c ---------------------------------------------------------------------
323 c First, calculate the limitation terms for NUT and Fe, and the
324 c Fe-limited Chl:C maximum. The light-saturated maximal photosynthesis
325 c rate term (Pc_m) is simply the product of a prescribed maximal
326 c photosynthesis rate (Pc_0), the Eppley temperature dependence, and a
327 c resource limitation term. The iron limitation term has a lower limit
328 c of Fe_lim_min and is scaled by (k_Fe2P + Fe2P_max) / Fe2P_max so that
329 c it approaches 1 as Fe approaches infinity. Thus, it is of comparable
330 c magnitude to the macro-nutrient limitation term.
331
332 c Macro-nutrient limitation
333 NO3_lim(i,j,k) = PTR_NO3(i,j,k)/(PTR_NO3(i,j,k)+k_NO3)
334
335 PO4_lim(i,j,k) = PTR_PO4(i,j,k)/(PTR_PO4(i,j,k)+k_PO4)
336
337 c Iron limitation
338
339 Fe_lim(i,j,k) = PTR_FE(i,j,k) / (PTR_FE(i,j,k)+k_Fe)
340
341 Fe_lim_diaz(i,j,k) = PTR_FE(i,j,k) / (PTR_FE(i,j,k)+k_Fe_diaz)
342
343 c ---------------------------------------------------------------------
344 c Diazotrophs are assumed to be more strongly temperature sensitive,
345 c given their observed restriction to relatively warm waters. Presumably
346 c this is because of the difficulty of achieving N2 fixation in an oxic
347 c environment. Thus, they have lower pc_0 and higher kappa_eppley.
348 c Taking the square root, to provide the geometric mean.
349
350 expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
351
352 c Light-saturated maximal photosynthesis rate
353
354 c Pc_m = Pc_0 * expkT(i,j,k)
355 c & * max(1. _d -8, NO3_lim(i,j,k) * PO4_lim(i,j,k)
356 c & * Fe_lim(i,j,k))**(1. / 3.)
357 c & * maskC(i,j,k,bi,bj)
358 c
359 c Pc_m_diaz = Pc_0_diaz
360 c & * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
361 c & * max(1. _d -8, PO4_lim(i,j,k)
362 c & * Fe_lim_diaz(i,j,k))**(1. / 2.)
363 c & * maskC(i,j,k,bi,bj)
364
365
366 #ifdef BLING_ADJOINT_SAFE_tmp_xxxxxxxxxxxxxxxxxx_needs_testing
367 th1 = tanh( (NO3_lim(i,j,k)-PO4_lim(i,j,k))*1. _d 6 )
368 nut_lim = ( 1. _d 0 - th1 ) * NO3_lim(i,j,k) * 0.5 _d 0
369 & + ( 1. _d 0 + th1 ) * PO4_lim(i,j,k) * 0.5 _d 0
370
371 th2 = tanh( (nut_lim-Fe_lim(i,j,k))*1. _d 6 )
372 tot_lim = ( 1. _d 0 - th2 ) * nut_lim * 0.5 _d 0
373 & + ( 1. _d 0 + th2 ) * Fe_lim(i,j,k) * 0.5 _d 0
374
375 th3 = tanh( (PO4_lim(i,j,k)-Fe_lim(i,j,k))*1. _d 6 )
376 diaz_lim = ( 1. _d 0 - th3 ) * PO4_lim(i,j,k) * 0.5 _d 0
377 & + ( 1. _d 0 + th3 ) * Fe_lim(i,j,k) * 0.5 _d 0
378
379
380 Pc_m = Pc_0 * expkT(i,j,k) * tot_lim
381 & * maskC(i,j,k,bi,bj)
382
383 Pc_m_diaz = Pc_0_diaz
384 & * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
385 & * diaz_lim * maskC(i,j,k,bi,bj)
386
387 #else
388
389 Pc_m = Pc_0 * expkT(i,j,k)
390 & * min(NO3_lim(i,j,k), PO4_lim(i,j,k), Fe_lim(i,j,k))
391 & * maskC(i,j,k,bi,bj)
392
393 Pc_m_diaz = Pc_0_diaz
394 & * exp(kappa_eppley_diaz * theta(i,j,k,bi,bj))
395 & * min(PO4_lim(i,j,k), Fe_lim_diaz(i,j,k))
396 & * maskC(i,j,k,bi,bj)
397
398 #endif
399
400
401 c ---------------------------------------------------------------------
402 c Fe limitation 1) reduces photosynthetic efficiency (alpha_Fe)
403 c and 2) reduces the maximum achievable Chl:C ratio (theta_Fe)
404 c below a prescribed, Fe-replete maximum value (theta_Fe_max),
405 c to approach a prescribed minimum Chl:C (theta_Fe_min) under extreme
406 c Fe-limitation.
407
408 theta_Fe_max = theta_Fe_max_lo+
409 & (theta_Fe_max_hi-theta_Fe_max_lo)*Fe_lim(i,j,k)
410
411 theta_Fe = theta_Fe_max/(1. _d 0 + alpha_photo*theta_Fe_max
412 & *irr_mem(i,j,k,bi,bj)/(epsln + 2. _d 0*Pc_m))
413
414 c ---------------------------------------------------------------------
415 c Nutrient-limited efficiency of algal photosystems, irrk, is calculated
416 c with the iron limitation term included as a multiplier of the
417 c theta_Fe_max to represent the importance of Fe in forming chlorophyll
418 c accessory antennae, which do not affect the Chl:C but still affect the
419 c phytoplankton ability to use light (eg Stzrepek & Harrison, Nature 2004).
420
421 irrk(i,j,k) = Pc_m/(epsln + alpha_photo*theta_Fe_max) +
422 & irr_mem(i,j,k,bi,bj)/2. _d 0
423
424 c Carbon-specific photosynthesis rate
425 mu(i,j,k) = Pc_m * ( 1. _d 0 - exp(-irr_eff(i,j,k)
426 & /(epsln + irrk(i,j,k))))
427
428 mu_diaz(i,j,k) = Pc_m_diaz * ( 1. _d 0 - exp(-irr_eff(i,j,k)
429 & /(epsln + irrk(i,j,k))))
430
431 ENDIF
432 ENDDO
433 ENDDO
434 ENDDO
435
436
437 C$TAF STORE P_sm = comlev1, key = ikey_dynamics, kind=isbyte
438 C$TAF STORE P_lg = comlev1, key = ikey_dynamics, kind=isbyte
439 C$TAF STORE P_diaz = comlev1, key = ikey_dynamics, kind=isbyte
440 Cxx needed?
441
442 c Instantaneous nutrient concentration in phyto biomass
443 c Separate loop so adjoint stuff above can be outside loop
444 c (fix for recomputations)
445
446 DO k=1,Nr
447 DO j=jmin,jmax
448 DO i=imin,imax
449
450 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
451
452 c expkT = exp(kappa_eppley * theta(i,j,k,bi,bj))
453
454 P_lg(i,j,k,bi,bj) = P_lg(i,j,k,bi,bj) +
455 & P_lg(i,j,k,bi,bj)*(mu(i,j,k) - lambda_0
456 & *expkT(i,j,k) *
457 & (P_lg(i,j,k,bi,bj)/pivotal)**(1. / 3.))
458 & * PTRACERS_dTLev(k)
459
460 P_sm(i,j,k,bi,bj) = P_sm(i,j,k,bi,bj) +
461 & P_sm(i,j,k,bi,bj)*(mu(i,j,k) - lambda_0
462 & *expkT(i,j,k) * (P_sm(i,j,k,bi,bj)/pivotal) )
463 & * PTRACERS_dTLev(k)
464
465 P_diaz(i,j,k,bi,bj) = P_diaz(i,j,k,bi,bj) +
466 & P_diaz(i,j,k,bi,bj)*(mu_diaz(i,j,k) - lambda_0
467 & *expkT(i,j,k) * (P_diaz(i,j,k,bi,bj)/pivotal) )
468 & * PTRACERS_dTLev(k)
469
470 ENDIF
471 ENDDO
472 ENDDO
473 ENDDO
474
475 C$TAF STORE P_sm = comlev1, key = ikey_dynamics, kind=isbyte
476 C$TAF STORE P_lg = comlev1, key = ikey_dynamics, kind=isbyte
477 C$TAF STORE P_diaz = comlev1, key = ikey_dynamics, kind=isbyte
478 cxx needed?
479
480
481 DO k=1,Nr
482 DO j=jmin,jmax
483 DO i=imin,imax
484
485 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
486
487 c use the diagnostic biomass to calculate the chl concentration
488 chl(i,j,k,bi,bj) = max(chl_min, CtoN * 12.01 * theta_Fe *
489 & (P_lg(i,j,k,bi,bj) + P_sm(i,j,k,bi,bj)
490 & + P_diaz(i,j,k,bi,bj)))
491
492 c stoichiometry
493 PtoN(i,j,k) = PtoN_min + (PtoN_max - PtoN_min) *
494 & PTR_PO4(i,j,k) / (k_PtoN + PTR_PO4(i,j,k))
495
496 FetoN(i,j,k) = FetoN_min + (FetoN_max - FetoN_min) *
497 & PTR_FE(i,j,k) / (k_FetoN + PTR_FE(i,j,k))
498
499 c Nutrient uptake
500 N_uptake(i,j,k) = mu(i,j,k)*(P_sm(i,j,k,bi,bj)
501 & + P_lg(i,j,k,bi,bj))
502
503 N_fix(i,j,k) = mu_diaz(i,j,k) * P_diaz(i,j,k,bi,bj)
504
505 P_uptake(i,j,k) = (N_uptake(i,j,k) +
506 & N_fix(i,j,k)) * PtoN(i,j,k)
507
508 Fe_uptake(i,j,k) = (N_uptake(i,j,k) +
509 & N_fix(i,j,k)) * FetoN(i,j,k)
510
511 c ---------------------------------------------------------------------
512 c Alkalinity is consumed through the production of CaCO3. Here, this is
513 c simply a linear function of the implied growth rate of small
514 c phytoplankton, which gave a reasonably good fit to the global
515 c observational synthesis of Dunne (2009). This is consistent
516 c with the findings of Jin et al. (GBC,2006).
517
518 CaCO3_uptake(i,j,k) = P_sm(i,j,k,bi,bj) * phi_sm *expkT(i,j,k)
519 & * mu(i,j,k) * CatoN
520
521 c ---------------------------------------------------------------------
522 c Partitioning between organic pools
523
524 c The uptake of nutrients is assumed to contribute to the growth of
525 c phytoplankton, which subsequently die and are consumed by heterotrophs.
526 c This can involve the transfer of nutrient elements between many
527 c organic pools, both particulate and dissolved, with complex histories.
528 c We take a simple approach here, partitioning the total uptake into two
529 c fractions - sinking and non-sinking - as a function of temperature,
530 c following Dunne et al. (2005).
531 c Then, the non-sinking fraction is further subdivided, such that the
532 c majority is recycled instantaneously to the inorganic nutrient pool,
533 c representing the fast turnover of labile dissolved organic matter via
534 c the microbial loop, and the remainder is converted to semi-labile
535 c dissolved organic matter. Iron and macro-nutrient are treated
536 c identically for the first step, but all iron is recycled
537 c instantaneously in the second step (i.e. there is no dissolved organic
538 c iron pool).
539
540 c sinking fraction: particulate organic matter
541
542 c expkT(i,j,k) = exp(kappa_eppley * theta(i,j,k,bi,bj))
543
544 frac_exp = (phi_sm + phi_lg * (mu(i,j,k)/
545 & (epsln + lambda_0*expkT(i,j,k)))**2.)/
546 & (1. + (mu(i,j,k)/(epsln + lambda_0*expkT(i,j,k)))**2.)*
547 & exp(kappa_remin * theta(i,j,k,bi,bj))
548
549 N_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
550 & (N_uptake(i,j,k) + N_fix(i,j,k))
551
552 P_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
553 & P_uptake(i,j,k)
554
555 Fe_spm(i,j,k) = frac_exp * (1.0 - phi_dvm) *
556 & Fe_uptake(i,j,k)
557
558 N_dvm(i,j,k) = frac_exp *
559 & (N_uptake(i,j,k) + N_fix(i,j,k)) - N_spm(i,j,k)
560
561 P_dvm(i,j,k) = frac_exp * P_uptake(i,j,k) -
562 & P_spm(i,j,k)
563
564 Fe_dvm(i,j,k) = frac_exp * Fe_uptake(i,j,k) -
565 & Fe_spm(i,j,k)
566
567 c the remainder is divided between instantaneously recycled and
568 c long-lived dissolved organic matter.
569
570 DON_prod(i,j,k) = phi_DOM*(N_uptake(i,j,k)
571 & + N_fix(i,j,k) - N_spm(i,j,k)
572 & - N_dvm(i,j,k))
573
574 DOP_prod(i,j,k) = phi_DOM*(P_uptake(i,j,k)
575 & - P_spm(i,j,k) - P_dvm(i,j,k))
576
577 N_recycle(i,j,k) = N_uptake(i,j,k) + N_fix(i,j,k)
578 & - N_spm(i,j,k) - DON_prod(i,j,k)
579 & - N_dvm(i,j,k)
580
581 P_recycle(i,j,k) = P_uptake(i,j,k)
582 & - P_spm(i,j,k) - DOP_prod(i,j,k)
583 & - P_dvm(i,j,k)
584
585 Fe_recycle(i,j,k) = Fe_uptake(i,j,k)
586 & - Fe_spm(i,j,k) - Fe_dvm(i,j,k)
587
588 ENDIF
589
590 ENDDO
591 ENDDO
592 ENDDO
593
594
595 c-----------------------------------------------------------
596 c remineralization of sinking organic matter
597 CALL BLING_REMIN(
598 I PTR_NO3, PTR_FE, PTR_O2, irr_inst,
599 I N_spm, P_spm, Fe_spm, CaCO3_uptake,
600 U N_reminp, P_reminp, Fe_reminsum,
601 U N_den_benthic, CACO3_diss,
602 I bi, bj, imin, imax, jmin, jmax,
603 I myIter, myTime, myThid)
604
605
606 c-----------------------------------------------------------
607 c remineralization from diel vertical migration
608 CALL BLING_DVM(
609 I N_dvm,P_dvm,Fe_dvm,
610 I PTR_O2, mld,
611 O N_remindvm, P_remindvm, Fe_remindvm,
612 I bi, bj, imin, imax, jmin, jmax,
613 I myIter, myTime, myThid)
614
615
616 c-----------------------------------------------------------
617 c sub grid scale sediments
618 #ifdef USE_SGS_SED
619 CALL BLING_SGS(
620 I xxx,
621 O xxx,
622 I bi, bj, imin, imax, jmin, jmax,
623 I myIter, myTime, myThid)#endif
624 #endif
625
626
627 c-----------------------------------------------------------
628 c
629
630 DO k=1,Nr
631 DO j=jmin,jmax
632 DO i=imin,imax
633
634 IF (hFacC(i,j,k,bi,bj) .gt. 0. _d 0) THEN
635
636
637 c Dissolved organic matter slow remineralization
638
639 #ifdef BLING_NO_NEG
640 DON_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DON
641 & *PTR_DON(i,j,k),0. _d 0)
642 DOP_remin(i,j,k) = MAX(maskC(i,j,k,bi,bj)*gamma_DOP
643 & *PTR_DOP(i,j,k),0. _d 0)
644 #else
645 DON_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DON
646 & *PTR_DON(i,j,k)
647 DOP_remin(i,j,k) = maskC(i,j,k,bi,bj)*gamma_DOP
648 & *PTR_DOP(i,j,k)
649 #endif
650
651
652 c Pelagic denitrification
653 c If anoxic
654 cxx IF (PTR_O2(i,j,k) .lt. 0. _d 0) THEN
655
656 IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
657 IF (PTR_NO3(i,j,k) .gt. oxic_min) THEN
658 N_den_pelag(i,j,k) = max(epsln, (NO3toN *
659 & ((1. _d 0 - phi_DOM) * (N_reminp(i,j,k)
660 & + N_remindvm(i,j,k)) + DON_remin(i,j,k) +
661 & N_recycle(i,j,k))) - N_den_benthic(i,j,k))
662 ENDIF
663 ENDIF
664
665 c Carbon flux diagnostic
666 POC_flux(i,j,k) = CtoN * N_spm(i,j,k)
667
668 NPP(i,j,k) = (N_uptake(i,j,k) + N_fix(i,j,k)) * CtoN
669
670 c oxygen production through photosynthesis
671 O2_prod(i,j,k) = O2toN * N_uptake(i,j,k)
672 & + (O2toN - 1.25 _d 0) * N_fix(i,j,k)
673
674
675
676 c-----------------------------------------------------------
677 C ADD TERMS
678
679 c Nutrients
680 c Sum of fast recycling, decay of sinking POM, and decay of DOM,
681 c less uptake, (less denitrification).
682
683 G_PO4(i,j,k) = -P_uptake(i,j,k) + P_recycle(i,j,k)
684 & + (1. _d 0 - phi_DOM) * (P_reminp(i,j,k)
685 & + P_remindvm(i,j,k)) + DOP_remin(i,j,k)
686
687 G_NO3(i,j,k) = -N_uptake(i,j,k)
688 IF (PTR_O2(i,j,k) .lt. oxic_min) THEN
689 c Anoxic
690 G_NO3(i,j,k) = G_NO3(i,j,k)
691 & - N_den_pelag(i,j,k) - N_den_benthic(i,j,k)
692 ELSE
693 c Oxic
694 G_NO3(i,j,k) = G_NO3(i,j,k)
695 & + N_recycle(i,j,k) + (1. _d 0 - phi_DOM) *
696 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
697 & + DON_remin(i,j,k)
698 ENDIF
699
700 cxxxx check
701 NCP(i,j,k) = (-G_NO3(i,j,k) + N_fix(i,j,k)) * CtoN
702
703 c Iron
704 c remineralization, sediments and adsorption are all bundled into
705 c Fe_reminsum
706
707 G_FE(i,j,k) = -Fe_uptake(i,j,k) + Fe_reminsum(i,j,k)
708 & + Fe_remindvm(i,j,k) + Fe_recycle(i,j,k)
709
710 c Dissolved Organic Matter
711 c A fraction of POM remineralization goes into dissolved pools.
712
713 G_DON(i,j,k) = DON_prod(i,j,k) + phi_DOM *
714 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
715 & - DON_remin(i,j,k)
716
717 G_DOP(i,j,k) = DOP_prod(i,j,k) + phi_DOM *
718 & (P_reminp(i,j,k) + P_remindvm(i,j,k))
719 & - DOP_remin(i,j,k)
720
721 c Oxygen:
722 c Assuming constant O2:N ratio in terms of oxidant required per mol of organic N.
723 c This implies a constant stoichiometry of C:N and H:N (where H is reduced, organic H).
724 c Because the N provided by N2 fixation is reduced from N2, rather than NO3-, the
725 c o2_2_n_fix is slightly less than the NO3- based ratio (by 1.25 mol O2/ mol N).
726 c Account for the organic matter respired through benthic denitrification by
727 c subtracting 5/4 times the benthic denitrification NO3 utilization rate from
728 c the overall oxygen consumption.
729
730 G_O2(i,j,k) = O2_prod(i,j,k)
731 c If oxic
732 IF (PTR_O2(i,j,k) .gt. oxic_min) THEN
733 G_O2(i,j,k) = G_O2(i,j,k)
734 & -O2toN * ((1. _d 0 - phi_DOM) *
735 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
736 & + DON_remin(i,j,k) + N_recycle(i,j,k))
737 c If anoxic but NO3 concentration is very low
738 c (generate negative O2; proxy for HS-).
739 ELSEIF (PTR_NO3(i,j,k) .lt. oxic_min) THEN
740 G_O2(i,j,k) = G_O2(i,j,k)
741 & -O2toN * ((1. _d 0 - phi_DOM) *
742 & (N_reminp(i,j,k) + N_remindvm(i,j,k))
743 & + DON_remin(i,j,k) + N_recycle(i,j,k))
744 & + N_den_benthic(i,j,k) * 1.25 _d 0
745 ENDIF
746
747 G_CaCO3(i,j,k) = CaCO3_diss(i,j,k) - CaCO3_uptake(i,j,k)
748 cxx sediments not accounted for
749
750 ENDIF
751
752 ENDDO
753 ENDDO
754 ENDDO
755
756
757 c ---------------------------------------------------------------------
758
759 #ifdef ALLOW_DIAGNOSTICS
760 IF ( useDiagnostics ) THEN
761
762 c 3d global variables
763 CALL DIAGNOSTICS_FILL(P_sm,'BLGPSM ',0,Nr,1,bi,bj,myThid)
764 CALL DIAGNOSTICS_FILL(P_lg,'BLGPLG ',0,Nr,1,bi,bj,myThid)
765 CALL DIAGNOSTICS_FILL(P_diaz,'BLGPDIA ',0,Nr,1,bi,bj,myThid)
766 CALL DIAGNOSTICS_FILL(chl,'BLGCHL ',0,Nr,1,bi,bj,myThid)
767 CALL DIAGNOSTICS_FILL(irr_mem,'BLGIMEM ',0,Nr,1,bi,bj,myThid)
768 c 3d local variables
769 CALL DIAGNOSTICS_FILL(irrk,'BLGIRRK ',0,Nr,2,bi,bj,myThid)
770 CALL DIAGNOSTICS_FILL(irr_eff,'BLGIEFF ',0,Nr,2,bi,bj,myThid)
771 CALL DIAGNOSTICS_FILL(Fe_lim,'BLGFELIM',0,Nr,2,bi,bj,myThid)
772 CALL DIAGNOSTICS_FILL(NO3_lim,'BLGNLIM ',0,Nr,2,bi,bj,myThid)
773 CALL DIAGNOSTICS_FILL(POC_flux,'BLGPOCF ',0,Nr,2,bi,bj,myThid)
774 CALL DIAGNOSTICS_FILL(NPP,'BLGNPP ',0,Nr,2,bi,bj,myThid)
775 CALL DIAGNOSTICS_FILL(NCP,'BLGNCP ',0,Nr,2,bi,bj,myThid)
776 c CALL DIAGNOSTICS_FILL(Fe_ads_inorg,'BLGFEAI',0,Nr,2,bi,bj,
777 c & myThid)
778 c CALL DIAGNOSTICS_FILL(Fe_dvm,'BLGFEDVM',0,Nr,2,bi,bj,myThid)
779 c CALL DIAGNOSTICS_FILL(Fe_sed,'BLGFESED',0,Nr,2,bi,bj,myThid)
780 CALL DIAGNOSTICS_FILL(Fe_spm,'BLGFESPM',0,Nr,2,bi,bj,myThid)
781 CALL DIAGNOSTICS_FILL(Fe_recycle,'BLGFEREC',0,Nr,2,bi,bj,
782 & myThid)
783 CALL DIAGNOSTICS_FILL(Fe_remindvm,'BLGFERD',0,Nr,2,bi,bj,
784 & myThid)
785 c CALL DIAGNOSTICS_FILL(Fe_reminp,'BLGFEREM',0,Nr,2,bi,bj,myThid)
786 CALL DIAGNOSTICS_FILL(Fe_reminsum,'BLGFEREM',0,Nr,2,bi,bj,
787 & myThid)
788 CALL DIAGNOSTICS_FILL(Fe_uptake,'BLGFEUP ',0,Nr,2,bi,bj,myThid)
789 CALL DIAGNOSTICS_FILL(N_den_benthic,'BLGNDENB',0,Nr,2,bi,bj,
790 & myThid)
791 CALL DIAGNOSTICS_FILL(N_den_pelag,'BLGNDENP',0,Nr,2,bi,bj,
792 & myThid)
793 CALL DIAGNOSTICS_FILL(N_dvm,'BLGNDVM ',0,Nr,2,bi,bj,myThid)
794 CALL DIAGNOSTICS_FILL(N_fix,'BLGNFIX ',0,Nr,2,bi,bj,myThid)
795 CALL DIAGNOSTICS_FILL(DON_prod,'BLGDONP ',0,Nr,2,bi,bj,myThid)
796 CALL DIAGNOSTICS_FILL(N_spm,'BLGNSPM ',0,Nr,2,bi,bj,myThid)
797 CALL DIAGNOSTICS_FILL(N_recycle,'BLGNREC ',0,Nr,2,bi,bj,myThid)
798 CALL DIAGNOSTICS_FILL(N_remindvm,'BLGNRD ',0,Nr,2,bi,bj,myThid)
799 CALL DIAGNOSTICS_FILL(N_reminp,'BLGNREM ',0,Nr,2,bi,bj,myThid)
800 CALL DIAGNOSTICS_FILL(N_uptake,'BLGNUP ',0,Nr,2,bi,bj,myThid)
801 CALL DIAGNOSTICS_FILL(P_dvm,'BLGPDVM ',0,Nr,2,bi,bj,myThid)
802 CALL DIAGNOSTICS_FILL(DOP_prod,'BLGDOPP ',0,Nr,2,bi,bj,myThid)
803 CALL DIAGNOSTICS_FILL(P_spm,'BLGPSPM ',0,Nr,2,bi,bj,myThid)
804 CALL DIAGNOSTICS_FILL(P_recycle,'BLGPREC ',0,Nr,2,bi,bj,myThid)
805 CALL DIAGNOSTICS_FILL(P_remindvm,'BLGPRD ',0,Nr,2,bi,bj,myThid)
806 CALL DIAGNOSTICS_FILL(P_reminp,'BLGPREM ',0,Nr,2,bi,bj,myThid)
807 CALL DIAGNOSTICS_FILL(P_uptake,'BLGPUP ',0,Nr,2,bi,bj,myThid)
808 c CALL DIAGNOSTICS_FILL(dvm,'BLGDVM ',0,Nr,2,bi,bj,myThid)
809 CALL DIAGNOSTICS_FILL(mu,'BLGMU ',0,Nr,2,bi,bj,myThid)
810 CALL DIAGNOSTICS_FILL(mu_diaz,'BLGMUDIA',0,Nr,2,bi,bj,myThid)
811 c 2d local variables
812 c CALL DIAGNOSTICS_FILL(Fe_burial,'BLGFEBUR',0,1,2,bi,bj,myThid)
813 c CALL DIAGNOSTICS_FILL(NO3_sed,'BLGNSED ',0,1,2,bi,bj,myThid)
814 c CALL DIAGNOSTICS_FILL(PO4_sed,'BLGPSED ',0,1,2,bi,bj,myThid)
815 c CALL DIAGNOSTICS_FILL(O2_sed,'BLGO2SED',0,1,2,bi,bj,myThid)
816 c these variables are currently 1d, could be 3d for diagnostics
817 c (or diag_fill could be called inside loop - which is faster?)
818 c CALL DIAGNOSTICS_FILL(frac_exp,'BLGFEXP ',0,Nr,2,bi,bj,myThid)
819 c CALL DIAGNOSTICS_FILL(irr_mix,'BLGIRRM ',0,Nr,2,bi,bj,myThid)
820 c CALL DIAGNOSTICS_FILL(irrk,'BLGIRRK ',0,Nr,2,bi,bj,myThid)
821 c CALL DIAGNOSTICS_FILL(kFe_eq_lig,'BLGPUP ',0,Nr,2,bi,bj,myThid)
822 c CALL DIAGNOSTICS_FILL(mu,'BLGMU ',0,Nr,2,bi,bj,myThid)
823 c CALL DIAGNOSTICS_FILL(mu_diaz,'BLGMUDIA',0,Nr,2,bi,bj,myThid)
824 c CALL DIAGNOSTICS_FILL(PtoN,'BLGP2N ',0,Nr,2,bi,bj,myThid)
825 c CALL DIAGNOSTICS_FILL(FetoN,'BLGFE2N ',0,Nr,2,bi,bj,myThid)
826 c CALL DIAGNOSTICS_FILL(Pc_m,'BLGPCM ',0,Nr,2,bi,bj,myThid)
827 c CALL DIAGNOSTICS_FILL(Pc_m_diaz,'BLGPCMD',0,Nr,2,bi,bj,myThid)
828 c CALL DIAGNOSTICS_FILL(theta_Fe,'BLGTHETA',0,Nr,2,bi,bj,myThid)
829 c CALL DIAGNOSTICS_FILL(theta_Fe_max,'BLGTHETM',0,Nr,2,bi,bj,myThid)
830 c CALL DIAGNOSTICS_FILL(wsink,'BLGWSINK',0,Nr,2,bi,bj,myThid)
831 c CALL DIAGNOSTICS_FILL(zremin,'BLGZREM ',0,Nr,2,bi,bj,myThid)
832 c CALL DIAGNOSTICS_FILL(z_dvm,'BLGZDVM ',0,Nr,2,bi,bj,myThid)
833
834 ENDIF
835 #endif /* ALLOW_DIAGNOSTICS */
836
837 #endif /* ALLOW_BLING */
838
839 RETURN
840
841 END

  ViewVC Help
Powered by ViewVC 1.1.22