/[MITgcm]/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F
ViewVC logotype

Contents of /MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F

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


Revision 1.15 - (show annotations) (download)
Thu May 1 16:19:32 2014 UTC (11 years, 2 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt65j_20150225, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt65_20140718, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt65e_20140929, ctrb_darwin2_ckpt64z_20140711, ctrb_darwin2_ckpt65l_20150504, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt64y_20140622, ctrb_darwin2_ckpt65d_20140915, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt64x_20140524, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt65g_20141120, ctrb_darwin2_ckpt65k_20150402, ctrb_darwin2_ckpt64w_20140502, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt65f_20141014, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt65i_20150123, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt65a_20140728, ctrb_darwin2_ckpt65b_20140812, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt65c_20140830, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt65h_20141217, ctrb_darwin2_ckpt66m_20171213, HEAD
Changes since 1.14: +35 -1 lines
add option FIX_ZOO_QUOTAS

1 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F,v 1.14 2013/06/12 17:48:20 jahn Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #include "PTRACERS_OPTIONS.h"
6 #include "DARWIN_OPTIONS.h"
7
8 #ifdef ALLOW_PTRACERS
9 #ifdef ALLOW_MONOD
10
11 c ====================================================================
12 c SUBROUTINE MONOD_PLANKTON
13 c 1. Local ecological interactions for models with many phytoplankton
14 c "functional groups"
15 c 2. Timestep plankton and nutrients locally
16 c 3. Includes explicit DOM and POM
17 c 4. Remineralization of detritus also determined in routine
18 c 5. Sinking particles and phytoplankton
19 c 6. NOT in this routine: iron chemistry
20 c
21 c Mick Follows, Scott Grant, Fall/Winter 2005
22 c Stephanie Dutkiewicz Spring/Summer 2006
23 c
24 c - add extra diagnostics, including R* (#define DAR_DIAG_RSTAR) - Stephanie, Spring 2007
25 c - add check for conservation (#define CHECK_CONS) - Stephanie, Spring 2007
26 c - improve grazing (#undef OLD_GRAZING) - Stephanie, Spring 2007
27 c - add diazotrophy (#define ALLOW_DIAZ) - Stephanie, Spring 2007
28 c - add mutation code (#define ALLOW_MUTANTS) - Jason Bragg, Spring/Summer 2007
29 c - new nitrogen limiting scheme (#undef OLD_NSCHEME) - Jason Bragg, Summer 2007
30 c - fix bug in diazotroph code - Stephanie, Fall 2007
31 c - add additional r* diagnostic for (no3+no2) - Stephanie, Winter 2007
32 c - add diversity diagnostics - Stephanie, Winter 2007
33 c - add geider chl:c ratio and growth rate dependence,
34 c though has no photo-inhibtion at this point - Stephanie, Spring 2008
35 c - add waveband dependence of light attenuation and absorption,
36 c NOTE: need to have geider turned on too - Anna Hickman, Summer 2008
37 c ====================================================================
38
39 c ANNA pass extra variables if WAVEBANDS
40 SUBROUTINE MONOD_PLANKTON(
41 U phyto,
42 I zooP, zooN, zooFe, zooSi,
43 O PP, Chl, Nfix, denit,
44 I PO4local, NO3local, FeTlocal, Silocal,
45 I NO2local, NH4local,
46 I DOPlocal, DONlocal, DOFelocal,
47 I POPlocal, PONlocal, POFelocal, PSilocal,
48 I phytoup, popuplocal, ponuplocal,
49 I pofeuplocal, psiuplocal,
50 I PARlocal,Tlocal, Slocal,
51 I pCO2local,
52 I freefelocal, inputFelocal,
53 I bottom, dzlocal,
54 O Rstarlocal, RNstarlocal,
55 #ifdef DAR_DIAG_GROW
56 O Growlocal, Growsqlocal,
57 #endif
58 #ifdef ALLOW_DIAZ
59 #ifdef DAR_DIAG_NFIXP
60 O NfixPlocal,
61 #endif
62 #endif
63 O dphytodt, dzooPdt, dzooNdt, dzooFedt,
64 O dzooSidt,
65 O dPO4dt, dNO3dt, dFeTdt, dSidt,
66 O dNH4dt, dNO2dt,
67 O dDOPdt, dDONdt, dDOFedt,
68 O dPOPdt, dPONdt, dPOFedt, dPSidt,
69 #ifdef ALLOW_CARBON
70 I DIClocal, DOClocal, POClocal, PIClocal,
71 I ALKlocal, O2local, ZooClocal,
72 I POCuplocal, PICuplocal,
73 O dDICdt, dDOCdt, dPOCdt, dPICdt,
74 O dALKdt, dO2dt, dZOOCdt,
75 #endif
76 #ifdef GEIDER
77 I phychl,
78 #ifdef DAR_DIAG_EK
79 O Ek, EkoverE,
80 #endif
81 #ifdef DAR_DIAG_PARW
82 O chl2c,
83 #endif
84 #ifdef DYNAMIC_CHL
85 O dphychl, Chlup,
86 #ifdef DAR_DIAG_EK
87 O acclim,
88 #endif
89 #endif
90 #ifdef ALLOW_CDOM
91 O dcdomdt ,
92 I cdomlocal,
93 #endif
94 #ifdef WAVEBANDS
95 I PARwlocal,
96 #ifdef DAR_DIAG_EK
97 O Ek_nl, EkoverE_nl,
98 #endif
99 #endif
100 #endif
101 #ifdef ALLOW_PAR_DAY
102 I PARdaylocal,
103 #endif
104 #ifdef DAR_DIAG_CHL
105 O ChlGeiderlocal, ChlDoneylocal,
106 O ChlCloernlocal,
107 #endif
108 I debug,
109 I runtim,
110 I MyThid)
111
112
113 implicit none
114 #include "MONOD_SIZE.h"
115 #include "MONOD.h"
116 #include "DARWIN_PARAMS.h"
117
118 c ANNA set wavebands params
119 #ifdef WAVEBANDS
120 #include "SPECTRAL_SIZE.h"
121 #include "WAVEBANDS_PARAMS.h"
122 #endif
123
124
125 C !INPUT PARAMETERS: ===================================================
126 C myThid :: thread number
127 INTEGER myThid
128 CEOP
129 c === GLOBAL VARIABLES =====================
130 c npmax = no of phyto functional groups
131 c nzmax = no of grazer species
132 c phyto = phytoplankton
133 c zoo = zooplankton
134 _RL phyto(npmax)
135 _RL zooP(nzmax)
136 _RL zooN(nzmax)
137 _RL zooFe(nzmax)
138 _RL zooSi(nzmax)
139 _RL PP
140 _RL Nfix
141 _RL denit
142 _RL Chl
143 _RL PO4local
144 _RL NO3local
145 _RL FeTlocal
146 _RL Silocal
147 _RL NO2local
148 _RL NH4local
149 _RL DOPlocal
150 _RL DONlocal
151 _RL DOFelocal
152 _RL POPlocal
153 _RL PONlocal
154 _RL POFelocal
155 _RL PSilocal
156 _RL phytoup(npmax)
157 _RL POPuplocal
158 _RL PONuplocal
159 _RL POFeuplocal
160 _RL PSiuplocal
161 _RL PARlocal
162 _RL Tlocal
163 _RL Slocal
164 _RL pCO2local
165 _RL freefelocal
166 _RL inputFelocal
167 _RL bottom
168 _RL dzlocal
169 _RL Rstarlocal(npmax)
170 _RL RNstarlocal(npmax)
171 #ifdef DAR_DIAG_GROW
172 _RL Growlocal(npmax)
173 _RL Growsqlocal(npmax)
174 #endif
175 #ifdef ALLOW_DIAZ
176 #ifdef DAR_DIAG_NFIXP
177 _RL NfixPlocal(npmax)
178 #endif
179 #endif
180 INTEGER debug
181 _RL dphytodt(npmax)
182 _RL dzooPdt(nzmax)
183 _RL dzooNdt(nzmax)
184 _RL dzooFedt(nzmax)
185 _RL dzooSidt(nzmax)
186 _RL dPO4dt
187 _RL dNO3dt
188 _RL dNO2dt
189 _RL dNH4dt
190 _RL dFeTdt
191 _RL dSidt
192 _RL dDOPdt
193 _RL dDONdt
194 _RL dDOFedt
195 _RL dPOPdt
196 _RL dPONdt
197 _RL dPOFedt
198 _RL dPSidt
199 #ifdef ALLOW_CARBON
200 _RL DIClocal
201 _RL DOClocal
202 _RL POClocal
203 _RL PIClocal
204 _RL ALKlocal
205 _RL O2local
206 _RL ZooClocal(nzmax)
207 _RL POCuplocal
208 _RL PICuplocal
209 _RL dDICdt
210 _RL dDOCdt
211 _RL dPOCdt
212 _RL dPICdt
213 _RL dALKdt
214 _RL dO2dt
215 _RL dZOOCdt(nzmax)
216 #endif
217 #ifdef GEIDER
218 _RL phychl(npmax)
219 _RL Ek(npmax)
220 _RL EkoverE(npmax)
221 _RL chl2c(npmax)
222 #ifdef DYNAMIC_CHL
223 _RL dphychl(npmax)
224 _RL Chlup(npmax)
225 _RL acclim(npmax)
226 #endif
227 #endif
228 #ifdef ALLOW_CDOM
229 _RL cdomlocal
230 _RL dcdomdt
231 #ifdef ALLOW_CARBON
232 _RL cdomclocal, dcdomcdt
233 #endif
234 #endif
235 #ifdef ALLOW_PAR_DAY
236 _RL PARdaylocal
237 #endif
238 #ifdef DAR_DIAG_CHL
239 _RL ChlGeiderlocal, ChlDoneylocal, ChlCloernlocal
240 #endif
241 _RL runtim
242
243 c ANNA Global variables for WAVEBANDS
244 c ANNA these variables are passed in/out of darwin_forcing.F
245 #ifdef WAVEBANDS
246 _RL PARwlocal(tlam) !PAR at midpoint of previous(in) and local(out) gridcell
247 _RL Ek_nl(npmax,tlam)
248 _RL EkoverE_nl(npmax,tlam)
249 #endif
250 c ANNA endif
251
252 c LOCAL VARIABLES
253 c -------------------------------------------------------------
254
255 c WORKING VARIABLES
256 c np = phytoplankton index
257 integer np
258 c nz = zooplankton index
259 integer nz
260
261 c variables for phytoplankton growth rate/nutrient limitation
262 c phytoplankton specific nutrient limitation term
263 _RL limit(npmax)
264 c phytoplankton light limitation term
265 _RL ilimit(npmax)
266 _RL pCO2limit(npmax)
267 _RL ngrow(npmax)
268 _RL grow(npmax)
269 _RL PspecificPO4(npmax)
270 _RL phytoTempFunction(npmax)
271 _RL mortPTempFunction
272 _RL dummy
273 _RL Ndummy
274 _RL Nsourcelimit(npmax)
275 _RL Nlimit(npmax)
276 _RL NO3limit(npmax)
277 _RL NO2limit(npmax)
278 _RL NH4limit(npmax)
279
280 c for check N limit scheme
281 _RL Ndiff
282 _RL NO3limcheck
283 _RL NO2limcheck
284 _RL Ndummy1
285 LOGICAL check_nlim
286
287 #ifndef OLD_NSCHEME
288 c [jbmodif] some new N terms
289 integer N2only
290 integer noNOdadv
291 integer NOreducost
292 _RL NO2zoNH4
293 _RL NOXzoNH4
294 #endif
295
296 c varible for mimumum phyto
297 _RL phytomin(npmax)
298
299 #ifdef OLD_GRAZE
300 c variables for zooplankton grazing rates
301 _RL zooTempFunction(nzmax)
302 _RL mortZTempFunction
303 _RL mortZ2TempFunction
304 _RL grazing_phyto(npmax)
305 _RL grazingP(nzmax)
306 _RL grazingN(nzmax)
307 _RL grazingFe(nzmax)
308 _RL grazingSi(nzmax)
309 #else
310 c variables for zooplankton grazing rates
311 _RL zooTempFunction(nzmax)
312 _RL mortZTempFunction
313 _RL mortZ2TempFunction
314 _RL allphyto(nzmax)
315 _RL denphyto(nzmax)
316 _RL grazphy(npmax,nzmax)
317 _RL sumgrazphy(npmax)
318 _RL sumgrazzoo(nzmax)
319 _RL sumgrazzooN(nzmax)
320 _RL sumgrazzooFe(nzmax)
321 _RL sumgrazzooSi(nzmax)
322 _RL sumgrazloss(nzmax)
323 _RL sumgrazlossN(nzmax)
324 _RL sumgrazlossFe(nzmax)
325 _RL sumgrazlossSi(nzmax)
326 #endif
327
328 #ifdef GEIDER
329 _RL alpha_I(npmax)
330 _RL pcarbon(npmax)
331 _RL pcm(npmax)
332 #ifdef DYNAMIC_CHL
333 _RL psinkchl(npmax)
334 _RL rhochl(npmax)
335 #endif
336 #endif
337
338 #ifdef ALLOW_CDOM
339 _RL cdomp_degrd, cdomn_degrd, cdomfe_degrd
340 _RL preminP_cdom, preminN_cdom, preminFe_cdom
341 #ifdef ALLOW_CARBON
342 _RL cdomc_degrd
343 _RL preminC_cdom
344 #endif
345 #endif
346
347 #ifdef DAR_DIAG_CHL
348 _RL tmppcm
349 _RL tmpchl2c
350 #endif
351 c variables for nutrient uptake
352 _RL consumpPO4
353 _RL consumpNO3
354 _RL consumpNO2
355 _RL consumpNH4
356 _RL consumpFeT
357 _RL consumpSi
358
359 c variables for reminerlaization of DOM and POM
360 _RL reminTempFunction
361 _RL DOPremin
362 _RL DONremin
363 _RL DOFeremin
364 _RL preminP
365 _RL preminN
366 _RL preminFe
367 _RL preminSi
368
369 c for sinking matter
370 _RL psinkP
371 _RL psinkN
372 _RL psinkFe
373 _RL psinkSi
374 _RL psinkphy(npmax)
375
376 #ifdef ALLOW_CARBON
377 _RL consumpDIC
378 _RL consumpDIC_PIC
379 _RL preminC
380 _RL DOCremin
381 _RL totphy_doc
382 _RL totzoo_doc
383 _RL totphy_poc
384 _RL totzoo_poc
385 _RL totphy_pic
386 _RL totzoo_pic
387 _RL psinkC
388 _RL psinkPIC
389 _RL disscPIC
390 #ifdef OLD_GRAZE
391 _RL grazingC(nzmax)
392 #else
393 c variables for zooplankton grazing rates
394 _RL sumgrazzooC(nzmax)
395 _RL sumgrazlossC(nzmax)
396 _RL sumgrazlossPIC(nzmax)
397 #endif
398
399 #endif
400
401 c variables for conversions from phyto and zoo to DOM and POM
402 _RL totphy_dop
403 _RL totphy_pop
404 _RL totphy_don
405 _RL totphy_pon
406 _RL totphy_dofe
407 _RL totphy_pofe
408 _RL totphy_dosi
409 _RL totphy_posi
410
411 _RL totzoo_dop
412 _RL totzoo_pop
413 _RL totzoo_don
414 _RL totzoo_pon
415 _RL totzoo_dofe
416 _RL totzoo_pofe
417 _RL totzoo_posi
418
419 #ifdef ALLOW_CDOM
420 _RL totphy_cdomp
421 _RL totphy_cdomn
422 _RL totphy_cdomfe
423 _RL totzoo_cdomp
424 _RL totzoo_cdomn
425 _RL totzoo_cdomfe
426 #ifdef ALLOW_CARBON
427 _RL totphy_cdomc
428 _RL totzoo_cdomc
429 #endif
430 #endif
431
432 _RL NO2prod
433 _RL NO3prod
434
435 _RL facpz
436
437 _RL kpar, kinh
438
439 _RL tmpr,tmpz, tmpgrow, tmp1, tmp2
440
441 integer ITEST
442
443 #ifdef PART_SCAV
444 _RL scav_part
445 _RL scav_poc
446 #endif
447
448
449 c ANNA local variables for WAVEBANDS
450 #ifdef WAVEBANDS
451 integer i,ilam
452 integer nl
453
454 c ANNA for interpolation
455 _RL cu_area
456 C _RL waves_diff
457 C _RL light_diff
458 C _RL alphaI_diff
459 C _RL squ_part
460 C _RL tri_part
461 C _RL seg_area
462
463 c ANNA inportant but local variables that can be fogotten
464 _RL PARwdn(tlam) !light at bottom of local gridcell
465 _RL attenwl(tlam) !attenuation (m-1)
466 _RL sumaphy_nl(tlam) !total phyto absorption at each wavelength
467 #endif
468 c ANNA endif
469
470 c.................................................................
471
472 #ifdef ALLOW_MUTANTS
473 c -m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-
474 c mutation variables [jbmodif]
475 INTEGER nsisone
476 INTEGER nsistwo
477 INTEGER nsisthree
478 INTEGER nsisfour
479 INTEGER npro
480 INTEGER taxind
481 _RL mutfor, mutback
482 _RL grow1
483 _RL grow2
484 _RL grow3
485 _RL grow4
486 #endif
487
488 INTEGER numtax
489 _RL oneyr,threeyr
490
491 #ifdef ALLOW_MUTANTS
492 c compile time options -- could maybe be moved to
493 c run time and set in data.gchem???
494 c QQQQQQQ
495 c Initialize sister taxon mutation scheme
496 c if numtax = 1, mutation is off
497 numtax = 4
498 c number of plankton types to assign for
499 c wild and mutants types
500 npro = 60
501 #else
502 numtax=1
503 #endif
504
505 oneyr = 86400.0 _d 0*360.0 _d 0
506 threeyr = oneyr*3. _d 0
507
508 c end mutation variables [jbmodif]
509 c -m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-
510
511 #ifndef OLD_NSCHEME
512 c [jbmodif] init new N terms
513 c if those not using NO3 has
514 c N limit with denominator with NO3 or not: 0=NO3 in denom; 1=NO2 only
515 N2only = 1
516 c ??
517 noNOdadv = 1
518 c energetic disadvantage of using NO2/No3: off=0, on=1
519 NOreducost =0
520 #endif
521
522 #ifdef GEIDER
523 do np=1,npmax
524 pcarbon(np) = 0. _d 0
525 pcm(np)=0. _d 0
526 chl2c(np)=0. _d 0
527 Ek(np)=0. _d 0
528 EkoverE(np)=0. _d 0
529 #ifdef DYNAMIC_CHL
530 acclim(np)=0. _d 0
531 psinkChl(np)=0. _d 0
532 #endif
533 #ifdef WAVEBANDS
534 do ilam=1,tlam
535 Ek_nl(np,ilam)=0. _d 0
536 EkoverE_nl(np,ilam)=0. _d 0
537 enddo
538 #endif
539 enddo
540 #endif
541
542
543 c set sum totals to zero
544 totphy_pop = 0. _d 0
545 totphy_dop = 0. _d 0
546 totphy_don = 0. _d 0
547 totphy_pon = 0. _d 0
548 totphy_dofe = 0. _d 0
549 totphy_pofe = 0. _d 0
550 totphy_posi = 0. _d 0
551
552 totzoo_dop = 0. _d 0
553 totzoo_pop = 0. _d 0
554 totzoo_don = 0. _d 0
555 totzoo_pon = 0. _d 0
556 totzoo_dofe = 0. _d 0
557 totzoo_pofe = 0. _d 0
558 totzoo_posi = 0. _d 0
559
560 consumpPO4 = 0.0 _d 0
561 consumpNO3 = 0.0 _d 0
562 consumpNO2 = 0.0 _d 0
563 consumpNH4 = 0.0 _d 0
564 consumpFeT = 0.0 _d 0
565 consumpSi = 0.0 _d 0
566
567 #ifdef ALLOW_CARBON
568 totphy_doc = 0. _d 0
569 totphy_poc = 0. _d 0
570 totphy_pic = 0. _d 0
571 totzoo_doc = 0. _d 0
572 totzoo_poc = 0. _d 0
573 totzoo_pic = 0. _d 0
574 consumpDIC = 0.0 _d 0
575 consumpDIC_PIC = 0.0 _d 0
576 #endif
577
578 c zeros for diagnostics
579 PP=0. _d 0
580 Nfix=0. _d 0
581 denit=0. _d 0
582 Chl=0. _d 0
583
584 c set up phtyoplankton array to be used for grazing and mortality
585 c set up other variable used more than once to zero
586 do np = 1, npmax
587 dummy = phyto(np)-phymin
588 phytomin(np)=max(dummy,0. _d 0)
589 NH4limit(np)=0. _d 0
590 NO2limit(np)=0. _d 0
591 NO3limit(np)=0. _d 0
592 #ifdef ALLOW_DIAZ
593 #ifdef DAR_DIAG_NFIXP
594 NfixPlocal(np)=0. _d 0
595 #endif
596 #endif
597 enddo
598
599
600 #ifdef ALLOW_MUTANTS
601 c SWD if parent population is zero (ie. negative) treat all mutants
602 c as zeros too
603 if(runtim .gt. threeyr) then
604 if(numtax .gt. 1)then
605 do np=1,npro
606 if(mod(np,numtax).eq. 1. _d 0)then
607 nsisone = np
608 nsistwo = np+1
609 nsisthree = np+2
610 nsisfour = np+3
611
612 if (phyto(nsisone).le.0. _d 0) then
613 if (numtax.gt.1) phyto(nsistwo)=0. _d 0
614 if (numtax.gt.2) phyto(nsisthree)=0. _d 0
615 if (numtax.gt.3) phyto(nsisfour)=0. _d 0
616 endif
617 endif
618 enddo
619 endif
620 endif
621 ccccccccccccccccccccccccccccccc
622 #endif
623
624
625 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
626 call MONOD_TEMPFUNC(Tlocal,phytoTempFunction,
627 & zooTempFunction, reminTempFunction,
628 & mortPTempFunction, mortZTempFunction,
629 & mortZ2TempFunction, myThid)
630 if (debug.eq.1) print*,'phytoTempFunction',
631 & phytoTempFunction, Tlocal
632 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
633
634 c ******************** GROWTH OF PHYTO ****************************
635 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
636 #ifndef GEIDER
637 c ANNA also if not wavebands
638 #ifndef WAVEBANDS
639 c Determine phytoplantkon light limitation: will affect growth rate
640 c using Platt-like equations with inhibition
641 do np = 1, npmax
642 if (PARlocal.gt.1. _d 0) then
643 kpar=ksatPAR(np)/10. _d 0;
644 kinh=kinhib(np)/1000. _d 0;
645 ilimit(np)=(1.0 _d 0 - EXP(-PARlocal*kpar))
646 & *(EXP(-PARlocal*kinh)) /
647 & ( kpar/(kpar+kinh)*EXP(kinh/kpar*LOG(kinh/(kpar+kinh))) )
648 ilimit(np)=min(ilimit(np),1. _d 0)
649 else
650 ilimit(np)=0. _d 0
651 endif
652 enddo
653 if (debug.eq.1) print*,'ilimit',ilimit, PARlocal
654 #endif
655 #endif
656 c ANNA endif
657
658 c pCO2 limit - default to non
659 do np=1,npmax
660 pCO2limit(np)=1. _d 0
661 c if (np.eq.6) then
662 c pCO2limit(np)=1. _d 0 + (pCO2local-400. _d -6)/600 _d -6
663 c pCO2limit(np)=max(pCO2limit(np),1. _d 0)
664 c pCO2limit(np)=min(pCO2limit(np),2. _d 0)
665 c endif
666 if (debug.eq.15) print*,'pco2limit',pCO2limit(np),pCO2local
667 enddo
668
669
670 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
671 c Determine phytoplankton nutrient limitation as mimimum of
672 c P,N,Si,Fe. However N can be utilized in several forms, so
673 c also determine which is used
674 do np=1, npmax
675 limit(np) = 1.0 _d 0
676 c P limitation
677 if (ksatPO4(np).gt.0. _d 0) then
678 dummy = PO4local/(PO4local+ksatPO4(np))
679 if (dummy .lt. limit(np)) limit(np) = dummy
680 endif
681 c Fe limitation
682 if (ksatFeT(np).gt.0. _d 0) then
683 dummy = FeTlocal/(FeTlocal+ksatFeT(np))
684 if (dummy .lt. limit(np))limit(np) = dummy
685 endif
686 c Si limiation
687 if (R_SiP(np) .ne. 0. _d 0.and.ksatSi(np).gt.0. _d 0) then
688 dummy = Silocal/(Silocal+ksatSi(np))
689 if (dummy .lt. limit(np))limit(np) = dummy
690 endif
691
692 c N limitation [jbmodif]
693 c nsource: genetic preference for {1:NH4&NO2 2:NH4 3:ALL Sources}
694 c Nsourcelimit marker for which nsource will be consumed {1:NO3 2:NO2 3:NH4}
695 c (Note: very different to way 1-D model does this)
696 if(diazotroph(np) .ne. 1.0 _d 0)then
697
698 c NH4, all nsource
699 if (ksatNH4(np).gt.0. _d 0) then
700 NH4limit(np) = NH4local/(NH4local+ksatNH4(np))
701 endif
702
703 #ifdef OLD_NSCHEME
704 if (ksatNO2(np).gt.0. _d 0) then
705 c NO2, if nsource is 1 or 3
706 NO2limit(np) = NO2local/(NO2local+ksatNO2(np))*
707 & EXP(-sig1*NH4local)
708 NO2limcheck = NO2local/(NO2local+ksatNO2(np))
709 endif
710 c NO3, if nsource is 3
711 if (ksatNO3(np).gt.0. _d 0) then
712 NO3limit(np) = NO3local/(NO3local+ksatNO3(np))*
713 & EXP(-sig2*NH4local - sig3*NO2local)
714 NO3limcheck = NO3local/(NO3local+ksatNO3(np))
715 endif
716 #else
717 c [jbmodif]
718 c NO2, if nsource is 1 or 3
719 if (ksatNO2(np).gt.0. _d 0 .and. nsource(np).ne.2) then
720 if (N2only.eq.1 .and. nsource(np).eq.1) then
721 c if (nsource(np).eq.1) then
722 NO2limit(np) = NO2local/(NO2local+ksatNO2(np))
723 & *EXP(-sig1*NH4local)
724 NO2limcheck = NO2local/(NO2local+ksatNO2(np))
725 else
726 if (ksatNO3(np).gt.0. _d 0) then
727 NO2limit(np)=NO2local/(NO3local+NO2local+ksatNO3(np))
728 & *EXP(-sig1*NH4local)
729 NO2limcheck=NO2local/(NO3local+NO2local+ksatNO3(np))
730 endif
731 endif
732 endif
733 c NO3, if nsource is 3
734 if (ksatNO3(np).gt.0. _d 0 .and. nsource(np).eq.3) then
735 NO3limit(np)=NO3local/(NO3local+NO2local+ksatNO3(np))
736 & *EXP(-sig1*NH4local)
737 NO3limcheck=NO3local/(NO3local+NO2local+ksatNO3(np))
738 endif
739
740 #endif
741
742 if (nsource(np).eq.2) then
743 NO2limit(np) = 0. _d 0
744 NO3limit(np) = 0. _d 0
745 NO2limcheck = 0. _d 0
746 NO3limcheck = 0. _d 0
747 endif
748 if (nsource(np).eq.1) then
749 NO3limit(np) = 0. _d 0
750 NO3limcheck = 0. _d 0
751 endif
752 if (nsource(np).eq.3) then
753 c don't do anything
754 endif
755
756 Ndummy = NO3limit(np)+NO2limit(np)+NH4limit(np)
757 c
758 c make sure no Nlim disadvantage;
759 c check that limit doesn't decrease at high NH4 levels
760 check_nlim=.FALSE.
761 if (check_nlim) then
762 Ndummy1=NO3limcheck+NO2limcheck
763 if (Ndummy.gt.0. _d 0.and.Ndummy.lt.Ndummy1) then
764 c print*,'QQ N limit WARNING',Ndummy, Ndummy1,
765 c & NO3local,NO2local,NH4local
766 Ndiff=Ndummy1-NH4limit(np)
767 NO2limit(np)=Ndiff *
768 & NO2limit(np)/(NO2limit(np)+NO3limit(np))
769 NO3limit(np)=Ndiff *
770 & NO3limit(np)/(NO2limit(np)+NO3limit(np))
771 Ndummy = NO3limit(np)+NO2limit(np)+NH4limit(np)
772 endif
773 endif
774
775 if (Ndummy.gt.1. _d 0) then
776 NO3limit(np) = NO3limit(np)/Ndummy
777 NO2limit(np) = NO2limit(np)/Ndummy
778 NH4limit(np) = NH4limit(np)/Ndummy
779 endif
780 Nlimit(np)=NO3limit(np)+NO2limit(np)+NH4limit(np)
781 if (Nlimit(np).gt.1.01 _d 0) then
782 print*,'QQ Nlimit', Nlimit(np), NO3limit(np),
783 & NO2limit(np), NH4limit(np)
784 endif
785 if (Nlimit(np).le.0. _d 0) then
786 c if (np.eq.1) then
787 c print*,'QQ Nlimit', Nlimit(np), NO3limit(np),
788 c & NO2limit(np), NH4limit(np)
789 c print*,'QQ limit',limit(np), np
790 c endif
791 Nlimit(np)=0. _d 0 !1 _d -10
792 endif
793
794 #ifdef OLD_NSCHEME
795 c lower growth for higher NO3 consumption at higher light
796 if (Nlimit(np).le.0. _d 0) then
797 ngrow(np)=1. _d 0
798 else
799 if (parlocal.gt.ilight) then
800 ngrow(np)=ngrowfac+(1. _d 0-ngrowfac)*
801 & (NH4limit(np)+NO2limit(np))/Nlimit(np)
802 else
803 ngrow(np)=1. _d 0
804 endif
805 ngrow(np)=min(ngrow(np),1. _d 0)
806 endif
807 #else
808 c disadvantage of oxidized inorganic N
809 c for now, ignore - a first attempt is included below
810 ngrow(np) = 1.0 _d 0
811
812 cc lower growth for higher NO3 consumption at higher light
813 c one possible way of counting cost of reducing NOX
814 if (NOreducost .eq. 1)then
815 if (Nlimit(np).le.0. _d 0) then
816 ngrow(np)=1. _d 0
817 else
818 ngrow(np)= (10. _d 0*4. _d 0 +2. _d 0) /
819 & (10. _d 0*4. _d 0 +2. _d 0*NH4limit(np)/Nlimit(np)
820 & +8. _d 0*NO2limit(np)/Nlimit(np)
821 & +10. _d 0*NO3limit(np)/Nlimit(np))
822 ngrow(np)=min(ngrow(np),1. _d 0)
823 endif
824 endif
825 c
826 c might consider other costs, too
827 c if (NOironcost .eq. 1)then
828 c
829 c endif
830 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
831 #endif
832
833 c Now Check Against General Nutrient Limiting Tendency
834 if (ksatNH4(np).gt.0. _d 0.or.ksatNO2(np).gt.0. _d 0
835 & .or.ksatNO3(np).gt.0. _d 0) then
836 if(Nlimit(np) .lt. limit(np)) limit(np) = Nlimit(np)
837 endif
838 else
839 ngrow(np)=1. _d 0
840 Nlimit(np)=1. _d 0
841 NO3limit(np)=0. _d 0
842 NO2limit(np)=0. _d 0
843 NH4limit(np)=0. _d 0
844 endif ! diaz
845 limit(np)=min(limit(np),1. _d 0)
846 enddo !np
847 if (debug.eq.1) print*,'nut limit',
848 & limit, PO4local, FeTlocal, Silocal
849 if (debug.eq.1) print*,'Nlimit',
850 & Nlimit
851 if (debug.eq.1) print*,'NH4limit',
852 & NH4limit, NH4local
853 if (debug.eq.1) print*,'NO2limit',
854 & NO2limit, NO2local
855 if (debug.eq.1) print*,'NO3limit',
856 & NO3limit, NO3local
857 if (debug.eq.1) print*,'ngrow',
858 & ngrow
859
860
861 #ifdef GEIDER
862
863 #ifdef WAVEBANDS
864 c ANNA if wavebands then uses spectral alphachl derived from spectral alpha * I
865 c so first get value for alphachl_nl * PARwlocal
866 c value will depend on matchup between spectra of alphachl_nl (ie. aphy_chl) and PARwlocal
867 c integrate alpha*PAR over wavebands
868 do np = 1,npmax
869 alpha_I(np) = 0 _d 0
870 do nl = 1,tlam
871 alpha_I(np) = alpha_I(np) + alphachl_nl(np,nl)*PARwlocal(nl)
872 end do
873 end do
874 c Geider growth (and chl2c) now depends on this (sinlge) value of alpha_chl * I
875
876 c alpha_mean now precomputed in darwin_init_vari
877 #else
878 c ANNA if not wavebands uses alphachl derived from mQyield * aphy_chl_ave
879 c for use with generic geider equation need to use alpha_I (ie. alphachl*PARlocal)
880 do np = 1, npmax
881 alpha_I(np)=alphachl(np)*PARlocal
882 enddo
883 c ANNA endif
884 #endif
885
886 do np = 1, npmax
887 pcm(np)=pcmax(np)*limit(np)*phytoTempFunction(np)*
888 & pCO2limit(np)
889 #ifdef DYNAMIC_CHL
890 if (phyto(np).gt. 0. _d 0) then
891 chl2c(np)=phychl(np)/(phyto(np)*R_PC(np))
892 else
893 chl2c(np)= 0. _d 0
894 endif
895 #endif
896 if (pcm(np).gt.0.d0) then
897 #ifndef DYNAMIC_CHL
898 c assumes balanced growth, eq A14 in Geider et al 1997
899 chl2c(np)=chl2cmax(np)/
900 & (1+(chl2cmax(np)*alpha_I(np))/
901 & (2*pcm(np)))
902 chl2c(np)=min(chl2c(np),chl2cmax(np))
903 chl2c(np)=max(chl2c(np),chl2cmin(np))
904 #endif
905 if (PARlocal.gt.1. _d -1) then
906 c Eq A1 in Geider et al 1997
907 pcarbon(np)=pcm(np)*( 1 -
908 & exp((-alpha_I(np)*chl2c(np))/(pcm(np))) )
909 #ifdef WAVEBANDS
910 Ek(np) = pcm(np)/(chl2c(np)*alpha_mean(np))
911 EkoverE(np) = Ek(np) / PARlocal
912 do nl=1, tlam
913 Ek_nl(np,nl)=pcm(np)/(chl2c(np)*alphachl_nl(np,nl))
914 EkoverE_nl(np,nl) = Ek_nl(np,nl) / PARwlocal(nl)
915 enddo
916 #else
917 Ek(np) = pcm(np)/(chl2c(np)*alphachl(np))
918 EkoverE(np) = Ek(np) / PARlocal
919 #endif
920 c for inhibition
921 if (inhibcoef_geid(np).gt.0. _d 0) then
922 if (PARlocal .ge. Ek(np)) then !photoinhibition begins
923 pcarbon(np) = pcarbon(np)*
924 & (EkoverE(np)*inhibcoef_geid(np))
925 endif
926 endif
927 c end inhib
928 if (pcarbon(np).lt. 0. _d 0)
929 & print*,'QQ ERROR pc=',np,pcarbon(np)
930 if (pcm(np).gt.0. _d 0) then
931 ilimit(np)=pcarbon(np)/pcm(np)
932 else
933 ilimit(np)= 0. _d 0
934 endif
935 else
936 ilimit(np)=0. _d 0
937 pcarbon(np)=0. _d 0
938 endif
939 else ! if pcm 0
940 pcm(np)=0. _d 0
941 #ifndef DYNAMIC_CHL
942 chl2c(np)=chl2cmin(np)
943 #endif
944 pcarbon(np)=0. _d 0
945 ilimit(np)=0. _d 0
946 endif
947 #ifdef DYNAMIC_CHL
948 c Chl:C acclimated to current conditions
949 c (eq A14 in Geider et al 1997)
950 acclim(np)=chl2cmax(np)/
951 & (1+(chl2cmax(np)*alpha_I(np))/
952 & (2*pcm(np)))
953 acclim(np)=min(acclim(np),chl2cmax(np))
954 c acclim(np)=max(acclim(np),chl2cmin(np))
955 #else
956 phychl(np)=phyto(np)*R_PC(np)*chl2c(np)
957 #endif
958 enddo
959 if (debug.eq.14) print*,'ilimit',ilimit, PARlocal
960 if (debug.eq.14) print*,'chl:c',chl2c
961 if (debug.eq.14) print*,'chl',phychl
962 #ifdef DYNAMIC_CHL
963 if (debug.eq.14) print*,'acclim',acclim
964 #endif
965 #endif /* GEIDER */
966
967 #ifdef DAR_DIAG_CHL
968 c diagnostic version of the above that does not feed back to growth
969 ChlGeiderlocal = 0. _d 0
970 do np = 1, npmax
971 tmppcm = mu(np)*limit(np)*phytoTempFunction(np)*
972 & pCO2limit(np)
973 if (tmppcm.gt.0.d0) then
974 tmpchl2c = Geider_chl2cmax(np)/
975 & (1+(Geider_chl2cmax(np)*Geider_alphachl(np)*PARdaylocal)/
976 & (2*tmppcm))
977 tmpchl2c = min(tmpchl2c, Geider_chl2cmax(np))
978 tmpchl2c = max(tmpchl2c, Geider_chl2cmin(np))
979 else
980 tmpchl2c = Geider_chl2cmin(np)
981 endif
982 ChlGeiderlocal = ChlGeiderlocal + phyto(np)*R_PC(np)*tmpchl2c
983 enddo
984 C Chl a la Doney
985 ChlDoneylocal = 0. _d 0
986 do np = 1, npmax
987 tmpchl2c = (Doney_Bmax - (Doney_Bmax-Doney_Bmin)*
988 & MIN(1. _d 0,PARdaylocal/Doney_PARstar))
989 & *limit(np)
990 ChlDoneylocal = ChlDoneylocal +
991 & tmpchl2c*R_PC(np)*phyto(np)
992 enddo
993 C Chl a la Cloern
994 ChlCloernlocal = 0. _d 0
995 do np = 1, npmax
996 tmpchl2c = Cloern_chl2cmin +
997 & Cloern_A*exp(Cloern_B*Tlocal)
998 & *exp(-Cloern_C*PARdaylocal)
999 & *limit(np)
1000 ChlCloernlocal = ChlCloernlocal +
1001 & tmpchl2c*R_PC(np)*phyto(np)
1002 enddo
1003 #endif /* DAR_DIAG_CHL */
1004
1005
1006 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1007 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1008 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1009 c ******************* END GROWTH PHYTO *******************************
1010
1011
1012 #ifdef OLD_GRAZE
1013 c------------------------------------------------------------------------
1014 c GRAZING sum contributions of all zooplankton
1015 do np=1,npmax
1016 grazing_phyto(np) = 0.0 _d 0
1017 do nz = 1, nzmax
1018 grazing_phyto(np) = grazing_phyto(np)
1019 & + graze(np,nz)*zooP(nz)*zooTempFunction(nz)
1020 enddo
1021 enddo
1022 if (debug.eq.2) print*,'grazing_phyto',grazing_phyto
1023 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1024 #else
1025 c------------------------------------------------------------------------
1026 c sum all palatability*phyto and find phyto specific grazing rate
1027 do nz=1,nzmax
1028 allphyto(nz)=0. _d 0
1029 do np=1,npmax
1030 allphyto(nz)=allphyto(nz)+palat(np,nz)*phyto(np)
1031 enddo
1032 if (allphyto(nz).le.0. _d 0) allphyto(nz)=phygrazmin
1033 #ifdef SER_GRAZ
1034 denphyto(nz)=0. _d 0
1035 do np=1,npmax
1036 denphyto(nz)=denphyto(nz)+
1037 & (palat(np,nz)*phyto(np)/allphyto(nz))*phyto(np)
1038 enddo
1039 if (denphyto(nz).le.0. _d 0) denphyto(nz)=phygrazmin
1040 #endif
1041 do np=1,npmax
1042 tmpz=max(0. _d 0,(allphyto(nz)-phygrazmin) )
1043 grazphy(np,nz)=grazemax(nz)*zooTempFunction(nz)*
1044 #ifdef SER_GRAZ
1045 c as in Vallina et al, 2011
1046 & (((palat(np,nz)*phyto(np)/allphyto(nz))*phyto(np))/
1047 & denphyto(nz)) *
1048 #else
1049 c as in Dutkiewicz et al, GBC, 2009
1050 & (palat(np,nz)*phyto(np)/allphyto(nz))*
1051 #endif
1052 & ( tmpz**hollexp/
1053 & (tmpz**hollexp+kgrazesat**hollexp) )
1054 enddo
1055 enddo
1056 if (debug.eq.2) print*,'allphyto',allphyto
1057 c if (debug.eq.2) print*,'grazephy',grazphy
1058 c sum over zoo for impact on phyto
1059 do np=1,npmax
1060 sumgrazphy(np)=0. _d 0
1061 do nz=1,nzmax
1062 sumgrazphy(np)=sumgrazphy(np)+
1063 & grazphy(np,nz)*zooP(nz)
1064 enddo
1065 enddo
1066 if (debug.eq.2) print*,'sumgrazephy',sumgrazphy
1067 c sum over phy for impact on zoo, and all remainder to go to POM
1068 do nz=1,nzmax
1069 sumgrazzoo(nz)=0. _d 0
1070 sumgrazzooN(nz)=0. _d 0
1071 sumgrazzooFe(nz)=0. _d 0
1072 sumgrazzooSi(nz)=0. _d 0
1073 sumgrazloss(nz)=0. _d 0
1074 sumgrazlossN(nz)=0. _d 0
1075 sumgrazlossFe(nz)=0. _d 0
1076 sumgrazlossSi(nz)=0. _d 0
1077 #ifdef ALLOW_CARBON
1078 sumgrazzooC(nz)=0. _d 0
1079 sumgrazlossC(nz)=0. _d 0
1080 sumgrazlossPIC(nz)=0. _d 0
1081 #endif
1082 do np=1,npmax
1083 #ifdef FIX_ZOO_QUOTAS
1084 sumgrazzoo(nz)=sumgrazzoo(nz)+
1085 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)
1086 sumgrazloss(nz)=sumgrazloss(nz)+
1087 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*zooP(nz)
1088 sumgrazzooN(nz)=sumgrazzooN(nz)+
1089 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_NP_zoo(nz)
1090 sumgrazlossN(nz)=sumgrazlossN(nz)+
1091 & (R_NP(np)-asseff(np,nz)*R_NP_zoo(nz))*
1092 & grazphy(np,nz)*zooP(nz)
1093 sumgrazzooFe(nz)=sumgrazzooFe(nz)+
1094 & asseff(np,nz)*grazphy(np,nz)*
1095 & zooP(nz)*R_FeP_zoo(nz)
1096 sumgrazlossFe(nz)=sumgrazlossFe(nz)+
1097 & (R_FeP(np)-asseff(np,nz)*R_FeP_zoo(nz))*
1098 & grazphy(np,nz)*zooP(nz)
1099 sumgrazzooSi(nz)=sumgrazzooSi(nz)+
1100 & asseff(np,nz)*grazphy(np,nz)*
1101 & zooP(nz)*R_SiP_zoo(nz)
1102 sumgrazlossSi(nz)=sumgrazlossSi(nz)+
1103 & (R_SiP(np)-asseff(np,nz)*R_SiP_zoo(nz))*
1104 & grazphy(np,nz)*zooP(nz)
1105 #ifdef ALLOW_CARBON
1106 sumgrazzooC(nz)=sumgrazzooC(nz)+
1107 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_PC_zoo(nz)
1108 sumgrazlossC(nz)=sumgrazlossC(nz)+
1109 & (R_PC(np)-asseff(np,nz)*R_PC_zoo(nz))*
1110 & grazphy(np,nz)*zooP(nz)
1111 sumgrazlossPIC(nz)=sumgrazlossPIC(nz)+
1112 & (1. _d 0)*grazphy(np,nz)*
1113 & zooP(nz)*R_PC(np)*R_PICPOC(np)
1114 #endif
1115 #else
1116 sumgrazzoo(nz)=sumgrazzoo(nz)+
1117 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)
1118 sumgrazloss(nz)=sumgrazloss(nz)+
1119 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*zooP(nz)
1120 sumgrazzooN(nz)=sumgrazzooN(nz)+
1121 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_NP(np)
1122 sumgrazlossN(nz)=sumgrazlossN(nz)+
1123 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1124 & zooP(nz)*R_NP(np)
1125 sumgrazzooFe(nz)=sumgrazzooFe(nz)+
1126 & asseff(np,nz)*grazphy(np,nz)*
1127 & zooP(nz)*R_FeP(np)
1128 sumgrazlossFe(nz)=sumgrazlossFe(nz)+
1129 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1130 & zooP(nz)*R_FeP(np)
1131 sumgrazzooSi(nz)=sumgrazzooSi(nz)+
1132 & asseff(np,nz)*grazphy(np,nz)*
1133 & zooP(nz)*R_SiP(np)
1134 sumgrazlossSi(nz)=sumgrazlossSi(nz)+
1135 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1136 & zooP(nz)*R_SiP(np)
1137 #ifdef ALLOW_CARBON
1138 sumgrazzooC(nz)=sumgrazzooC(nz)+
1139 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_PC(np)
1140 sumgrazlossC(nz)=sumgrazlossC(nz)+
1141 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1142 & zooP(nz)*R_PC(np)
1143 sumgrazlossPIC(nz)=sumgrazlossPIC(nz)+
1144 & (1. _d 0)*grazphy(np,nz)*
1145 & zooP(nz)*R_PC(np)*R_PICPOC(np)
1146 #endif
1147 #endif
1148 enddo
1149 enddo
1150 if (debug.eq.2) print*,'sumgrazzoo',sumgrazzoo
1151 if (debug.eq.2) print*,'sumgrazloss',sumgrazloss
1152 if (debug.eq.2) print*,'sumgrazzooN',sumgrazzooN
1153 if (debug.eq.2) print*,'sumgrazlossN',sumgrazlossN
1154 if (debug.eq.2) print*,'sumgrazzooFe',sumgrazzooFe
1155 if (debug.eq.2) print*,'sumgrazlossFe',sumgrazlossFe
1156 if (debug.eq.2) print*,'sumgrazzooSi',sumgrazzooSi
1157 if (debug.eq.2) print*,'sumgrazlossSi',sumgrazlossSi
1158 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1159 #endif
1160
1161 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1162 c accumulate particulate and dissolved detritus
1163 do np=1, npmax
1164 totphy_pop=totphy_pop+
1165 & ExportFracP(np)*mortphy(np)*
1166 & mortPTempFunction*phytomin(np)
1167 totphy_dop=totphy_dop+
1168 & (1. _d 0-ExportFracP(np))*mortphy(np)*
1169 & mortPTempFunction*phytomin(np)
1170 totphy_pon=totphy_pon+ R_NP(np)*
1171 & ExportFracP(np)*mortphy(np)*
1172 & mortPTempFunction*phytomin(np)
1173 totphy_don=totphy_don+ R_NP(np)*
1174 & (1. _d 0-ExportFracP(np))*mortphy(np)*
1175 & mortPTempFunction*phytomin(np)
1176 totphy_pofe=totphy_pofe+ R_FeP(np)*
1177 & ExportFracP(np)*mortphy(np)*
1178 & mortPTempFunction*phytomin(np)
1179 totphy_dofe=totphy_dofe+ R_FeP(np)*
1180 & (1. _d 0-ExportFracP(np))*mortphy(np)*
1181 & mortPTempFunction*phytomin(np)
1182 totphy_posi=totphy_posi+ R_SiP(np)*
1183 & mortphy(np)*
1184 & mortPTempFunction*phytomin(np)
1185 #ifdef ALLOW_CARBON
1186 totphy_poc=totphy_poc+ R_PC(np)*
1187 & ExportFracP(np)*mortphy(np)*
1188 & mortPTempFunction*phytomin(np)
1189 totphy_doc=totphy_doc+ R_PC(np)*
1190 & (1. _d 0-ExportFracP(np))*mortphy(np)*
1191 & mortPTempFunction*phytomin(np)
1192 totphy_pic=totphy_pic+ R_PC(np)*R_PICPOC(np)*
1193 & mortphy(np)*
1194 & mortPTempFunction*phytomin(np)
1195 #endif
1196 enddo
1197 #ifdef ALLOW_CDOM
1198 totphy_cdomp=(fraccdom)*totphy_dop
1199 totphy_dop=totphy_dop-totphy_cdomp
1200 totphy_cdomn=rnp_cdom*totphy_cdomp
1201 totphy_don=totphy_don-totphy_cdomn
1202 totphy_cdomfe=rfep_cdom*totphy_cdomp
1203 totphy_dofe=totphy_dofe-totphy_cdomfe
1204 #ifdef ALLOW_CARBON
1205 totphy_cdomc=rcp_cdom*totphy_cdomp
1206 totphy_doc=totphy_doc-totphy_cdomc
1207 #endif
1208 #endif
1209 if (debug.eq.3) print*,'tot_phy_pop',totphy_pop
1210 if (debug.eq.3) print*,'tot_phy_dop',totphy_dop
1211 if (debug.eq.3) print*,'tot_phy_pon',totphy_pon
1212 if (debug.eq.3) print*,'tot_phy_don',totphy_don
1213 if (debug.eq.3) print*,'tot_phy_pofe',totphy_pofe
1214 if (debug.eq.3) print*,'tot_phy_dofe',totphy_dofe
1215 if (debug.eq.3) print*,'tot_phy_posi',totphy_posi
1216
1217
1218 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1219
1220
1221 #ifdef OLD_GRAZE
1222 c ****************** ZOO GRAZING RATE ****************************
1223 c determine zooplankton grazing rates
1224 do nz = 1, nzmax
1225 c grazing: sum contribution from all phytoplankton
1226 grazingP(nz) = 0.0 _d 0
1227 grazingN(nz) = 0.0 _d 0
1228 grazingFe(nz) = 0.0 _d 0
1229 grazingSi(nz) = 0.0 _d 0
1230 #ifdef ALLOW_CARBON
1231 grazingC(nz) = 0.0 _d 0
1232 #endif
1233 do np = 1, npmax
1234 facpz = (phytomin(np)/(phytomin(np) + kgrazesat))
1235 & *zooTempFunction(nz)
1236 grazingP(nz) = grazingP(nz) +
1237 & graze(np,nz)*facpz
1238 grazingN(nz) = grazingN(nz) +
1239 & graze(np,nz)*R_NP(np)*facpz
1240 grazingFe(nz) = grazingFe(nz) +
1241 & graze(np,nz)*R_FeP(np)*facpz
1242 grazingSi(nz) = grazingSi(nz) +
1243 & graze(np,nz)*R_SiP(np)*facpz
1244 #ifdef ALLOW_CARBON
1245 grazingC(nz) = grazingC(nz) +
1246 & graze(np,nz)*R_PC(np)*facpz
1247 #endif
1248 enddo
1249 enddo
1250 if (debug.eq.4) print*,'grazingP', grazingP
1251 if (debug.eq.4) print*,'grazingN', grazingN
1252 if (debug.eq.4) print*,'grazingFe', grazingFe
1253 if (debug.eq.4) print*,'grazingSi', grazingSi
1254 c ************* END ZOO GRAZING *********************************
1255 #endif
1256 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1257 c accumulate particulate and dissolved detritus
1258 do nz=1, nzmax
1259 totzoo_pop=totzoo_pop+
1260 & ExportFracZ(nz)*( mortzoo(nz)*
1261 & mortZTempFunction*zooP(nz)
1262 & + mortzoo2(nz)*
1263 & mortZ2TempFunction*zooP(nz)**2 )
1264 totzoo_dop=totzoo_dop+
1265 & (1. _d 0-ExportFracZ(nz))*(
1266 & mortzoo(nz)*
1267 & mortZTempFunction*zooP(nz)+
1268 & mortzoo2(nz)*
1269 & mortZ2TempFunction*zooP(nz)**2 )
1270 totzoo_pon=totzoo_pon+
1271 & ExportFracZ(nz)*( mortzoo(nz)*
1272 & mortZTempFunction*zooN(nz)
1273 & + mortzoo2(nz)*
1274 & mortZ2TempFunction*zooN(nz)*zooP(nz) )
1275 totzoo_don=totzoo_don+
1276 & (1. _d 0-ExportFracZ(nz))*(
1277 & mortzoo(nz)*
1278 & mortZTempFunction*zooN(nz)+
1279 & mortzoo2(nz)*
1280 & mortZ2TempFunction*zooN(nz)*zooP(nz) )
1281 totzoo_pofe=totzoo_pofe+
1282 & ExportFracZ(nz)*( mortzoo(nz)*
1283 & mortZTempFunction*zooFe(nz)
1284 & + mortzoo2(nz)*
1285 & mortZ2TempFunction*zooFe(nz)*zooP(nz) )
1286 totzoo_dofe=totzoo_dofe+
1287 & (1. _d 0-ExportFracZ(nz))*(
1288 & mortzoo(nz)*
1289 & mortZTempFunction*zooFe(nz) +
1290 & mortzoo2(nz)*
1291 & mortZ2TempFunction*zooFe(nz)*zooP(nz) )
1292 totzoo_posi=totzoo_posi+
1293 & ( mortzoo(nz)*
1294 & mortZTempFunction*zooSi(nz)+
1295 & mortzoo2(nz)*
1296 & mortZ2TempFunction*zooSi(nz)*zooP(nz) )
1297 #ifdef ALLOW_CARBON
1298 totzoo_poc=totzoo_poc+
1299 & ExportFracZ(nz)*( mortzoo(nz)*
1300 & mortZTempFunction*zooClocal(nz)
1301 & + mortzoo2(nz)*
1302 & mortZ2TempFunction*zooClocal(nz)*zooP(nz) )
1303 totzoo_doc=totzoo_doc+
1304 & (1. _d 0-ExportFracZ(nz))*( mortzoo(nz)*
1305 & mortZTempFunction*zooClocal(nz)
1306 & + mortzoo2(nz)*
1307 & mortZ2TempFunction*zooClocal(nz)*zooP(nz) )
1308 #endif
1309 enddo
1310
1311 #ifndef OLD_GRAZE
1312 do nz=1, nzmax
1313 totzoo_pop=totzoo_pop+
1314 & ExportFracGraz(nz)*sumgrazloss(nz)
1315 totzoo_dop=totzoo_dop+
1316 & (1. _d 0-ExportFracGraz(nz))*sumgrazloss(nz)
1317 totzoo_pon=totzoo_pon+
1318 & ExportFracGraz(nz)*sumgrazlossN(nz)
1319 totzoo_don=totzoo_don+
1320 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossN(nz)
1321 totzoo_pofe=totzoo_pofe+
1322 & ExportFracGraz(nz)*sumgrazlossFe(nz)
1323 totzoo_dofe=totzoo_dofe+
1324 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossFe(nz)
1325 totzoo_posi=totzoo_posi+
1326 & sumgrazlossSi(nz)
1327 #ifdef ALLOW_CARBON
1328 totzoo_poc=totzoo_poc+
1329 & ExportFracGraz(nz)*sumgrazlossC(nz)
1330 totzoo_doc=totzoo_doc+
1331 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossC(nz)
1332 totzoo_pic=totzoo_pic+
1333 & sumgrazlossPIC(nz)
1334 #endif
1335 enddo
1336 #endif
1337
1338 #ifdef ALLOW_CDOM
1339 totzoo_cdomp=(fraccdom)*totzoo_dop
1340 totzoo_dop=totzoo_dop-totzoo_cdomp
1341 totzoo_cdomn=rnp_cdom*totzoo_cdomp
1342 totzoo_don=totzoo_don-totzoo_cdomn
1343 totzoo_cdomfe=rfep_cdom*totzoo_cdomp
1344 totzoo_dofe=totzoo_dofe-totzoo_cdomfe
1345 #ifdef ALLOW_CARBON
1346 totzoo_cdomc=rcp_cdom*totzoo_cdomp
1347 totzoo_doc=totzoo_doc-totzoo_cdomc
1348 #endif
1349 #endif
1350 if (debug.eq.5) print*,'totzoo_pop',totzoo_pop
1351 if (debug.eq.5) print*,'totzoo_dop',totzoo_dop
1352 if (debug.eq.5) print*,'totzoo_pon',totzoo_pon
1353 if (debug.eq.5) print*,'totzoo_don',totzoo_don
1354 if (debug.eq.5) print*,'totzoo_pofe',totzoo_pofe
1355 if (debug.eq.5) print*,'totzoo_dofe',totzoo_dofe
1356 if (debug.eq.5) print*,'totzoo_posi',totzoo_posi
1357 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1358
1359 c ********************* NUTRIENT UPTAKE *******************************
1360 c determine nutrient uptake
1361 c consumption - sum of phytoplankton contributions
1362 do np = 1, npmax
1363 c phospate uptake by each phytoplankton
1364 #ifndef GEIDER
1365 grow(np)=ngrow(np)*mu(np)*limit(np)*ilimit(np)*
1366 & phytoTempFunction(np)*pCO2limit(np)
1367 #endif
1368 #ifdef GEIDER
1369 grow(np)=ngrow(np)*pcarbon(np)
1370 if (debug.eq.1) print*,'grow', grow(np), pcarbon(np)
1371 if (debug.eq.14) print*,'grow', grow(np), pcarbon(np)
1372 #ifdef DYNAMIC_CHL
1373 c geider 97 for dChl/dt (source part) Eq. 3
1374 if (acclim(np).gt. 0. _d 0.and.
1375 & alpha_I(np).gt. 0. _d 0) then
1376 rhochl(np)=chl2cmax(np) *
1377 & (grow(np)/(alpha_I(np)*acclim(np)) )
1378 else
1379 rhochl(np)= 0. _d 0
1380 endif
1381 if (debug.eq.14) print*,'rhochl',rhochl(np)
1382 #endif
1383 #endif
1384 PspecificPO4(np) = grow(np)*phyto(np)
1385 c write(6,*)'np =',np, ' PspecificPO4 ='
1386 c & ,PspecificPO4(np)
1387 consumpPO4 = consumpPO4 + PspecificPO4(np)
1388 consumpFeT = consumpFeT + PspecificPO4(np)*R_FeP(np)
1389 consumpSi = consumpSi + PspecificPO4(np)*R_SiP(np)
1390 cswd should have O2prod as function of np?
1391 c New Way of doing Nitrogen Consumption .......................
1392 if(diazotroph(np) .ne. 1.0 _d 0)then
1393 if (Nlimit(np).le.0. _d 0) then
1394 consumpNO3 = consumpNO3
1395 consumpNO2 = consumpNO2
1396 consumpNH4 = consumpNH4
1397 else
1398 consumpNO3 = consumpNO3 +
1399 & NO3limit(np)/Nlimit(np)*PspecificPO4(np)*R_NP(np)
1400 consumpNO2 = consumpNO2 +
1401 & NO2limit(np)/Nlimit(np)* PspecificPO4(np)*R_NP(np)
1402 consumpNH4 = consumpNH4 +
1403 & NH4limit(np)/Nlimit(np)*PspecificPO4(np)*R_NP(np)
1404 endif
1405 else
1406 consumpNO3 = consumpNO3
1407 consumpNO2 = consumpNO2
1408 consumpNH4 = consumpNH4
1409 Nfix=Nfix+PspecificPO4(np)*R_NP(np)
1410 #ifdef ALLOW_DIAZ
1411 #ifdef DAR_DIAG_NFIXP
1412 NfixPlocal(np)=PspecificPO4(np)*R_NP(np)
1413 #endif
1414 #endif
1415 endif
1416 #ifdef ALLOW_CARBON
1417 consumpDIC = consumpDIC + PspecificPO4(np)*R_PC(np)
1418 consumpDIC_PIC = consumpDIC_PIC +
1419 & PspecificPO4(np)*R_PC(np)*R_PICPOC(np)
1420 #endif
1421 enddo
1422 if (debug.eq.7) print*,'local', parlocal,tlocal,po4local,
1423 & no3local, no2local,nh4local,fetlocal,silocal
1424 if (debug.eq.7) print*,'grow',grow
1425 if (debug.eq.6) print*,'pspecificpo4', PspecificPO4
1426 if (debug.eq.6) print*,'consumpPO4', consumpPO4
1427 if (debug.eq.6) print*,'consumpFeT', consumpFeT
1428 if (debug.eq.6) print*,'consumpSi ', consumpsi
1429 if (debug.eq.6) print*,'consumpNO3', consumpNO3
1430 if (debug.eq.6) print*,'consumpNO2', consumpNO2
1431 if (debug.eq.6) print*,'consumpNH4', consumpNH4
1432 c ****************** END NUTRIENT UPTAKE ****************************
1433
1434 c sinking phytoplankton and POM
1435 if(bottom .eq. 1.0 _d 0)then
1436 psinkP = (wp_sink*POPuplocal)/(dzlocal)
1437 psinkN = (wn_sink*PONuplocal)/(dzlocal)
1438 psinkFe = (wfe_sink*POFeuplocal)/(dzlocal)
1439 psinkSi = (wsi_sink*PSiuplocal)/(dzlocal)
1440 do np=1,npmax
1441 psinkPhy(np) =
1442 & (wsink(np)*Phytoup(np))/(dzlocal)
1443 enddo
1444 #ifdef DYNAMIC_CHL
1445 do np=1,npmax
1446 psinkChl(np) =
1447 & (wsink(np)*Chlup(np))/(dzlocal)
1448 enddo
1449 #endif
1450 #ifdef ALLOW_CARBON
1451 psinkC = (wc_sink*POCuplocal)/(dzlocal)
1452 psinkPIC = (wpic_sink*PICuplocal)/(dzlocal)
1453 #endif
1454 else
1455 psinkP = (wp_sink*(POPuplocal-POPlocal))/(dzlocal)
1456 psinkN = (wn_sink*(PONuplocal-PONlocal))/(dzlocal)
1457 psinkFe = (wfe_sink*(POFeuplocal-POFelocal))/(dzlocal)
1458 psinkSi = (wsi_sink*(PSiuplocal-PSilocal))/(dzlocal)
1459 do np=1,npmax
1460 psinkPhy(np) =
1461 & (wsink(np))*(Phytoup(np)-Phyto(np))/(dzlocal)
1462 enddo
1463 #ifdef DYNAMIC_CHL
1464 do np=1,npmax
1465 psinkChl(np) =
1466 & (wsink(np))*(Chlup(np)-phychl(np))/(dzlocal)
1467 enddo
1468 #endif
1469 #ifdef ALLOW_CARBON
1470 psinkC = (wc_sink*(POCuplocal-POClocal))/(dzlocal)
1471 psinkPIC = (wpic_sink*(PICuplocal-PIClocal))/(dzlocal)
1472 #endif
1473 endif
1474
1475 c DOM remineralization rates
1476 DOPremin = reminTempFunction * Kdop * DOPlocal
1477 DONremin = reminTempFunction * Kdon * DONlocal
1478 DOFeremin = reminTempFunction * KdoFe * DOFelocal
1479
1480 c remineralization of sinking particulate
1481 preminP = reminTempFunction * Kpremin_P*POPlocal
1482 preminN = reminTempFunction * Kpremin_N*PONlocal
1483 preminFe = reminTempFunction * Kpremin_Fe*POFelocal
1484 preminSi = reminTempFunction * Kpremin_Si*PSilocal
1485 #ifdef ALLOW_CDOM
1486 preminP_cdom = fraccdom*preminP
1487 preminP=preminP-preminP_cdom
1488 preminN_cdom = rnp_cdom*preminP_cdom
1489 preminN=preminN-preminN_cdom
1490 preminFe_cdom = rfep_cdom*preminP_cdom
1491 preminFe=preminFe-preminFe_cdom
1492 #endif
1493
1494 #ifdef ALLOW_CARBON
1495 DOCremin = reminTempFunction * Kdoc * DOClocal
1496 preminC = reminTempFunction * Kpremin_C*POClocal
1497 #ifdef ALLOW_CDOM
1498 preminC_cdom = rcp_cdom*preminP_cdom
1499 preminC=preminC-preminC_cdom
1500 #endif
1501 c dissolution
1502 disscPIC = Kdissc*PIClocal
1503 #endif
1504
1505 #ifdef ALLOW_CDOM
1506 c degradation of CDOM - high when bleached by light
1507 cdomp_degrd = reminTempFunction * cdomlocal
1508 & *(cdomdegrd+cdombleach*min(PARlocal/PARcdom,1. _d 0) )
1509 cdomn_degrd = rnp_cdom * cdomp_degrd
1510 cdomfe_degrd= rfep_cdom * cdomp_degrd
1511 #ifdef ALLOW_CARBON
1512 cdomc_degrd = rcp_cdom * cdomp_degrd
1513 #endif
1514 #endif
1515
1516 #ifdef ALLOW_DENIT
1517 if (O2local.lt.O2crit) then
1518 if (NO3local.lt.no3crit) then
1519 c no remineralization for N, P, Fe (not Si?)
1520 DOPremin = 0. _d 0
1521 DONremin = 0. _d 0
1522 DOFeremin = 0. _d 0
1523 preminP = 0. _d 0
1524 preminN = 0. _d 0
1525 preminFe = 0. _d 0
1526 #ifdef ALLOW_CDOM
1527 preminP_cdom = 0. _d 0
1528 preminN_cdom = 0. _d 0
1529 preminFe_cdom = 0. _d 0
1530 #endif
1531 #ifdef ALLOW_CARBON
1532 DOCremin = 0. _d 0
1533 preminC = 0. _d 0
1534 #ifdef ALLOW_CDOM
1535 preminC_cdom = 0. _d 0
1536 #endif
1537 #endif
1538
1539 #ifdef ALLOW_CDOM
1540 c degradation of CDOM - high when bleached by light
1541 cdomp_degrd = reminTempFunction * cdomlocal
1542 & *(cdombleach*min(PARlocal/PARcdom,1. _d 0) )
1543 cdomn_degrd = rnp_cdom * cdomp_degrd
1544 cdomfe_degrd= rfep_cdom * cdomp_degrd
1545 #ifdef ALLOW_CARBON
1546 cdomc_degrd = rcp_cdom * cdomp_degrd
1547 #endif
1548 #endif
1549 endif
1550 endif
1551 #endif
1552 c end denit caveats
1553
1554 c chemistry
1555 c NH4 -> NO2 -> NO3 by bacterial action
1556 NO2prod = knita*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1557 & *NH4local
1558 NO3prod = knitb*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1559 & *NO2local
1560 c NO2prod = knita*NH4local
1561 c NO3prod = knitb*NO2local
1562 c
1563 #ifdef PART_SCAV
1564 scav_poc=POPlocal/1.1321 _d -4
1565 c scav rate
1566 scav_part=scav_rat*scav_inter*(scav_poc**scav_exp)
1567 #endif
1568 c -------------------------------------------------------------------
1569 c calculate tendency terms (and some diagnostics)
1570 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1571 c phytoplankton
1572 do np=1,npmax
1573 dphytodt(np) = PspecificPO4(np)
1574 #ifdef OLD_GRAZE
1575 & - grazing_phyto(np)*
1576 & (phytomin(np)/(phytomin(np) + kgrazesat))
1577 #else
1578 & - sumgrazphy(np)
1579 #endif
1580 & - mortphy(np)*
1581 & mortPTempFunction*phytomin(np)
1582 & + psinkphy(np)
1583 #ifdef GEIDER
1584 #ifdef DYNAMIC_CHL
1585 dphychl(np) = acclim(np)*PspecificPO4(np)*R_PC(np)
1586 c dphychl(np) = rhochl(np)*PspecificPO4(np)*R_PC(np)
1587 & + acclimtimescl *
1588 & (acclim(np)-chl2c(np))*phyto(np)*R_PC(np)
1589 & +(
1590 #ifdef OLD_GRAZE
1591 & - grazing_phyto(np)*
1592 & (phytomin(np)/(phytomin(np) + kgrazesat))
1593 #else
1594 & - sumgrazphy(np)
1595 #endif
1596 & - mortphy(np)*
1597 & mortPTempFunction*phytomin(np))
1598 & *chl2c(np)*R_PC(np)
1599 & + psinkChl(np)
1600 #endif
1601 Chl=Chl + phychl(np)
1602 #endif
1603 c %% diagnostics
1604 PP = PP + PspecificPO4(np)
1605 c%%%
1606 #ifdef OLD_GRAZE
1607 tmpr=grazing_phyto(np)*
1608 & (phytomin(np)/(phytomin(np) + kgrazesat))
1609 & + mortphy(np)*
1610 & mortPTempFunction*phytomin(np)
1611 & - psinkphy(np)
1612 #else
1613 tmpr=sumgrazphy(np)
1614 & + mortphy(np)*
1615 & mortPTempFunction*phytomin(np)
1616 & - psinkphy(np)
1617 #endif
1618 #ifdef DAR_DIAG_RSTAR
1619 #ifndef GEIDER
1620 tmpgrow=ngrow(np)*mu(np)*ilimit(np)*
1621 & phytoTempFunction(np)*pCO2limit(np)
1622 #endif
1623 #ifdef GEIDER
1624 tmpgrow=grow(np)/limit(np)
1625 #endif
1626 tmp1=tmpgrow*phyto(np)-tmpr
1627 tmp2=tmpgrow*phyto(np)*(exp(-sig1*nh4local)+NH4limit(np))
1628 & -tmpr
1629 if (tmp1.ne.0. _d 0) then
1630 Rstarlocal(np)=ksatPO4(np)*tmpr/tmp1
1631 else
1632 Rstarlocal(np)=-9999. _d 0
1633 endif
1634 if (tmp2.ne.0. _d 0) then
1635 RNstarlocal(np)=ksatNO3(np)*
1636 & (tmpr-tmpgrow*NH4limit(np)*phyto(np))/tmp2
1637 else
1638 RNstarlocal(np)=-9999. _d 0
1639 endif
1640 #endif
1641 #ifdef DAR_DIAG_GROW
1642 c include temp, light, nutrients
1643 c Growlocal(np)=grow(np)
1644 c include temp and light, but not nutrients
1645 Growlocal(np)=ngrow(np)*mu(np)*ilimit(np)*
1646 & phytoTempFunction(np)
1647 c include temp, but not nutrients or light
1648 c Growlocal(np)=ngrow(np)*mu(np)*
1649 c & phytoTempFunction(np)
1650 Growsqlocal(np)=Growlocal(np)**2
1651 #endif
1652 enddo
1653 c end np loop
1654 if (debug.eq.10) print*,'dphytodt',dphytodt
1655 c
1656 #ifdef OLD_GRAZE
1657 c zooplankton growth by grazing
1658 do nz=1,nzmax
1659 c zoo in P currency
1660 dzooPdt(nz) = grazingP(nz)*zooP(nz)
1661 C zooplankton stoichiometry varies according to food source
1662 dzooNdt(nz) = grazingN(nz)*zooP(nz)
1663 dzooFedt(nz) = grazingFe(nz)*zooP(nz)
1664 dzooSidt(nz) = grazingSi(nz)*zooP(nz)
1665 enddo
1666 #else
1667 do nz=1,nzmax
1668 c zoo in P currency
1669 dzooPdt(nz) = sumgrazzoo(nz)
1670 C zooplankton stoichiometry varies according to food source
1671 dzooNdt(nz) = sumgrazzooN(nz)
1672 dzooFedt(nz) = sumgrazzooFe(nz)
1673 dzooSidt(nz) = sumgrazzooSi(nz)
1674 enddo
1675 #endif
1676 if (debug.eq.10) print*,'dZooPdt',dZooPdt
1677
1678 c zooplankton mortality
1679 do nz=1,nzmax
1680 c zoo in P currency
1681 dzooPdt(nz) = dzooPdt(nz)
1682 & - mortzoo(nz)*
1683 & mortZTempFunction*zooP(nz)
1684 & - mortzoo2(nz)*
1685 & mortZ2TempFunction*zooP(nz)**2
1686 c zooplankton in other currencies
1687 C zooplankton stoichiometry varies according to food source
1688 dzooNdt(nz) = dzooNdt(nz)
1689 & - mortzoo(nz)*
1690 & mortZTempFunction*zooN(nz)
1691 & - mortzoo2(nz)*
1692 & mortZ2TempFunction*zooN(nz)*zooP(nz)
1693 dzooFedt(nz) = dzooFedt(nz)
1694 & - mortzoo(nz)*
1695 & mortZTempFunction*zooFe(nz)
1696 & - mortzoo2(nz)*
1697 & mortZ2TempFunction*zooFe(nz)*zooP(nz)
1698 dzooSidt(nz) = dzooSidt(nz)
1699 & - mortzoo(nz)*
1700 & mortZTempFunction*zooSi(nz)
1701 & - mortzoo2(nz)*
1702 & mortZ2TempFunction*zooSi(nz)*zooP(nz)
1703 enddo
1704
1705
1706 c sum contributions to inorganic nutrient tendencies
1707 dPO4dt = - consumpPO4 +
1708 #ifdef ALLOW_CDOM
1709 & DOPremin
1710 #else
1711 & preminP + DOPremin
1712 #endif
1713 dNH4dt = - consumpNH4 +
1714 #ifdef ALLOW_CDOM
1715 & DONremin
1716 #else
1717 & preminN + DONremin
1718 #endif
1719 & - NO2prod
1720 dNO2dt = - consumpNO2
1721 & + NO2prod - NO3prod
1722 dNO3dt = - consumpNO3
1723 & + NO3prod
1724 c-ONLYNO3 dNO3dt = - consumpNO3 + preminN + DONremin
1725 #ifdef ALLOW_DENIT
1726 if (O2local.le.O2crit) then
1727 denit = denit_np
1728 #ifdef ALLOW_CDOM
1729 & *(DOPremin)
1730 #else
1731 & *(preminP + DOPremin)
1732 #endif
1733 dNO3dt = dNO3dt - (denit_no3/denit_np)*denit
1734 dNH4dt = dNH4dt -
1735 #ifdef ALLOW_CDOM
1736 & (DONremin)
1737 #else
1738 & (preminN + DONremin)
1739 #endif
1740 endif
1741 #endif
1742 dFeTdt = - consumpFeT +
1743 #ifdef ALLOW_CDOM
1744 & DOFeremin
1745 #else
1746 & preminFe + DOFeremin
1747 #endif
1748 #ifdef PART_SCAV
1749 & - scav_part*freefelocal +
1750 #else
1751 & - scav*freefelocal +
1752 #endif
1753 & alpfe*inputFelocal/dzlocal
1754 dSidt = - consumpSi + preminSi
1755
1756 c tendency of dissolved organic pool
1757 dDOPdt = totphy_dop + totzoo_dop - DOPremin
1758 #ifdef ALLOW_CDOM
1759 & +preminP + cdomp_degrd
1760 #endif
1761 dDONdt = totphy_don + totzoo_don - DONremin
1762 #ifdef ALLOW_CDOM
1763 & +preminN + cdomn_degrd
1764 #endif
1765 dDOFedt = totphy_dofe + totzoo_dofe - DOFeremin
1766 #ifdef ALLOW_CDOM
1767 & +preminFe + cdomfe_degrd
1768 #endif
1769
1770 c tendency of particulate detritus pools
1771 dpopdt = totphy_pop + totzoo_pop - preminP + psinkP
1772 #ifdef ALLOW_CDOM
1773 & -preminP_cdom
1774 #endif
1775 dpondt = totphy_pon + totzoo_pon - preminN + psinkN
1776 #ifdef ALLOW_CDOM
1777 & -preminN_cdom
1778 #endif
1779 dpofedt = totphy_pofe + totzoo_pofe - preminFe + psinkFe
1780 #ifdef ALLOW_CDOM
1781 & -preminFe_cdom
1782 #endif
1783 dpSidt = totphy_posi + totzoo_posi - preminSi + psinkSi
1784 #ifdef ALLOW_CARBON
1785 dDICdt = - consumpDIC - consumpDIC_PIC
1786 #ifdef ALLOW_CDOM
1787 & + DOCremin
1788 #else
1789 & + preminC + DOCremin
1790 #endif
1791 & + disscPIC
1792 dDOCdt = totphy_doc + totzoo_doc - DOCremin
1793 #ifdef ALLOW_CDOM
1794 & +preminC + cdomc_degrd
1795 #endif
1796 dPOCdt = totphy_poc + totzoo_poc - preminC + psinkC
1797 #ifdef ALLOW_CDOM
1798 & -preminC_cdom
1799 #endif
1800 dPICdt = totphy_pic + totzoo_pic - disscPIC + psinkPIC
1801 dALKdt = - dNO3dt - 2.d0 * (consumpDIC_PIC - disscPIC)
1802 c should be = O2prod - preminP - DOPremin?
1803 c OLD WAY
1804 c dO2dt = - R_OP*dPO4dt
1805 c production of O2 by photosynthesis
1806 dO2dt = R_OP*consumpPO4
1807 c loss of O2 by remineralization
1808 if (O2local.gt.O2crit) then
1809 dO2dt = dO2dt - R_OP
1810 #ifdef ALLOW_CDOM
1811 & *(DOPremin)
1812 #else
1813 & *(preminP + DOPremin)
1814 #endif
1815 endif
1816 #ifdef OLD_GRAZE
1817 do nz=1,nzmax
1818 dzooCdt(nz) = grazingC(nz)*zooClocal(nz)
1819 & - mortzoo(nz)*
1820 & mortZTempFunction*zooClocal(nz)
1821 & - mortzoo2(nz)*
1822 & mortZ2TempFunction*zooClocal(nz)*zooP(nz)
1823 enddo
1824 #else
1825 do nz=1,nzmax
1826 dzooCdt(nz) = sumgrazzooc(nz)
1827 & - mortzoo(nz)*
1828 & mortZTempFunction*zooClocal(nz)
1829 & - mortzoo2(nz)*
1830 & mortZ2TempFunction*zooClocal(nz)*zooP(nz)
1831 enddo
1832 #endif
1833
1834 #endif /* ALLOW_CARBON */
1835
1836 #ifdef ALLOW_CDOM
1837 dcdomdt = totphy_cdomp + totzoo_cdomp + preminP_cdom
1838 & -cdomp_degrd
1839 #endif
1840
1841 if (debug.eq.10) print*,'dDOPdt', dDOPdt
1842 if (debug.eq.10) print*,'dpopdt',dpopdt
1843 if (debug.eq.10) print*,'dDONdt',dDONdt
1844 if (debug.eq.10) print*,'dpondt',dpondt
1845 c
1846 c -------------------------------------------------------------------
1847 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1848 c --------------------------------------------------------------------------
1849
1850 c -m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-
1851 c Mutation - apply mutation to tendencies [jbmodif]
1852
1853 #ifdef ALLOW_MUTANTS
1854 c apply to all sisters when first sister is encountered
1855 if(runtim .gt. threeyr) then
1856 mutfor=1 _d -8
1857 mutback=1 _d -12
1858 if(numtax .gt. 1)then
1859 do np=1,npro
1860 if(mod(np,numtax).eq. 1. _d 0)then
1861 nsisone = np
1862 nsistwo = np+1
1863 nsisthree = np+2
1864 nsisfour = np+3
1865
1866 grow1 = PspecificPO4(nsisone)
1867 grow2 = PspecificPO4(nsistwo)
1868
1869 if(numtax.eq.2)grow3 = 0.0 _d 0
1870 if(numtax.eq.2)grow4 = 0.0 _d 0
1871
1872 if(numtax.eq.3)grow4 = 0.0 _d 0
1873 if(numtax.ge.3)grow3 = PspecificPO4(nsisthree)
1874
1875 if(numtax.eq.4)grow4 = PspecificPO4(nsisfour)
1876
1877
1878
1879 dphytodt(nsisone) = dphytodt(nsisone)
1880 & - grow1 *1.4427 _d 0*mutfor
1881 & - grow1 *1.4427 _d 0*mutfor
1882 & - grow1 *1.4427 _d 0*mutfor
1883 & + grow2 *1.4427 _d 0*mutback
1884 & + grow3 *1.4427 _d 0*mutback
1885 & + grow4 *1.4427 _d 0*mutback
1886
1887 dphytodt(nsistwo) = dphytodt(nsistwo)
1888 & - grow2 *1.4427 _d 0*mutback
1889 & + grow1 *1.4427 _d 0*mutfor
1890
1891 if(numtax .ge. 3)then
1892 dphytodt(nsisthree) = dphytodt(nsisthree)
1893 & - grow3 *1.4427 _d 0*mutback
1894 & + grow1 *1.4427 _d 0*mutfor
1895 endif
1896
1897 if(numtax .eq. 4)then
1898 dphytodt(nsisfour) = dphytodt(nsisfour)
1899 & - grow4 *1.4427 _d 0*mutback
1900 & + grow1 *1.4427 _d 0*mutfor
1901 c QQQQQQQQQQ FIX FOR NIT RUNS ONLY!!!
1902 if (phyto(nsisfour).eq.0. _d 0) then
1903 if (phyto(nsistwo).eq.0. _d 0) then
1904 if (dphytodt(nsistwo).gt.dphytodt(nsisfour)) then
1905 dphytodt(nsisfour)=dphytodt(nsistwo)
1906 endif
1907 endif
1908 if (phyto(nsisthree).eq.0. _d 0) then
1909 if (dphytodt(nsisthree).gt.dphytodt(nsisfour)) then
1910 dphytodt(nsisfour)=dphytodt(nsisthree)
1911 endif
1912 endif
1913 endif
1914 c QQQQQQQQQQQQQ
1915 endif
1916
1917 c QQQQQQQQQQQQTEST
1918 if (debug.eq.11) then
1919 if (PARlocal.gt.1. _d 0) then
1920 if (dphytodt(nsistwo).gt.dphytodt(nsisfour).and.
1921 & dphytodt(nsisfour).gt.0. _d 0) then
1922 print*,'QQQQ nsistwo>nsisfour',nsistwo,nsisfour,
1923 & dphytodt(nsistwo), dphytodt(nsisfour),
1924 & phyto(nsistwo), phyto(nsisfour),
1925 & phyto(nsisone)
1926 endif
1927 if (dphytodt(nsisthree).gt.dphytodt(nsisfour).and.
1928 & dphytodt(nsisfour).gt.0. _d 0) then
1929 print*,'QQQQ nsisthree>nsisfour',nsisthree,nsisfour,
1930 & dphytodt(nsisthree), dphytodt(nsisfour),
1931 & phyto(nsisthree), phyto(nsisfour),
1932 & phyto(nsisone)
1933 endif
1934 if (dphytodt(nsisfour).gt.dphytodt(nsisone).and.
1935 & dphytodt(nsisone).gt.0. _d 0) then
1936 print*,' BIG QQQQ nsisfour>nsisone',nsisone,nsisfour,
1937 & dphytodt(nsisfour), dphytodt(nsisone),
1938 & phyto(nsisfour), phyto(nsisone)
1939 endif
1940 endif
1941 endif
1942 c QQQQQQQQQTEST
1943 endif
1944 enddo
1945 endif
1946 endif
1947
1948 c mutation is finished
1949 c -m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-m-
1950 #endif
1951
1952
1953
1954 RETURN
1955 END
1956 #endif /*MONOD*/
1957 #endif /*ALLOW_PTRACERS*/
1958 c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22