/[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.9 - (show annotations) (download)
Thu Jul 26 18:01:22 2012 UTC (13 years ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt63r_20120817, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt64_20121012
Changes since 1.8: +40 -2 lines
o bug fix, so denitrification (and remineralization) will not happen when
  NO3 below a critical value (NO3crit)

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

  ViewVC Help
Powered by ViewVC 1.1.22