/[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.6 - (show annotations) (download)
Thu May 31 21:08:25 2012 UTC (13 years, 2 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63n_20120604
Changes since 1.5: +152 -7 lines
o add CDOM-like tracer (#define ALLOW_CDOM)

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

  ViewVC Help
Powered by ViewVC 1.1.22