/[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.3 - (show annotations) (download)
Wed May 4 23:22:25 2016 UTC (9 years, 3 months ago) by mmazloff
Branch: MAIN
Changes since 1.2: +17 -20 lines
mods for adjoint stability

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

  ViewVC Help
Powered by ViewVC 1.1.22