/[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.8 - (show annotations) (download)
Fri Jun 29 20:41:59 2012 UTC (13 years, 1 month ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63p_20120707
Changes since 1.7: +22 -5 lines
o add potential for CO2 limitation to growth rate -- commented out for now
  (need to check monod_forcing and monod_plankton to make it work)

1 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F,v 1.7 2012/06/28 20:36:11 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 c chemistry
1460 c NH4 -> NO2 -> NO3 by bacterial action
1461 NO2prod = knita*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1462 & *NH4local
1463 NO3prod = knitb*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1464 & *NO2local
1465 c NO2prod = knita*NH4local
1466 c NO3prod = knitb*NO2local
1467 c
1468 #ifdef PART_SCAV
1469 scav_poc=POPlocal/1.1321 _d -4
1470 c scav rate
1471 scav_part=scav_rat*scav_inter*(scav_poc**scav_exp)
1472 #endif
1473 c -------------------------------------------------------------------
1474 c calculate tendency terms (and some diagnostics)
1475 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1476 c phytoplankton
1477 do np=1,npmax
1478 dphytodt(np) = PspecificPO4(np)
1479 #ifdef OLD_GRAZE
1480 & - grazing_phyto(np)*
1481 & (phytomin(np)/(phytomin(np) + kgrazesat))
1482 #else
1483 & - sumgrazphy(np)
1484 #endif
1485 & - mortphy(np)*
1486 & mortPTempFunction*phytomin(np)
1487 & + psinkphy(np)
1488 #ifdef GEIDER
1489 #ifdef DYNAMIC_CHL
1490 dphychl(np) = acclim(np)*PspecificPO4(np)*R_PC(np)
1491 c dphychl(np) = rhochl(np)*PspecificPO4(np)*R_PC(np)
1492 & + acclimtimescl *
1493 & (acclim(np)-chl2c(np))*phyto(np)*R_PC(np)
1494 & +(
1495 #ifdef OLD_GRAZE
1496 & - grazing_phyto(np)*
1497 & (phytomin(np)/(phytomin(np) + kgrazesat))
1498 #else
1499 & - sumgrazphy(np)
1500 #endif
1501 & - mortphy(np)*
1502 & mortPTempFunction*phytomin(np))
1503 & *chl2c(np)*R_PC(np)
1504 & + psinkChl(np)
1505 #endif
1506 Chl=Chl + phychl(np)
1507 #endif
1508 c %% diagnostics
1509 PP = PP + PspecificPO4(np)
1510 c%%%
1511 #ifdef OLD_GRAZE
1512 tmpr=grazing_phyto(np)*
1513 & (phytomin(np)/(phytomin(np) + kgrazesat))
1514 & + mortphy(np)*
1515 & mortPTempFunction*phytomin(np)
1516 & - psinkphy(np)
1517 #else
1518 tmpr=sumgrazphy(np)
1519 & + mortphy(np)*
1520 & mortPTempFunction*phytomin(np)
1521 & - psinkphy(np)
1522 #endif
1523 #ifdef DAR_DIAG_RSTAR
1524 #ifndef GEIDER
1525 tmpgrow=ngrow(np)*mu(np)*ilimit(np)*
1526 & phytoTempFunction(np)*pCO2limit(np)
1527 #endif
1528 #ifdef GEIDER
1529 tmpgrow=grow(np)/limit(np)
1530 #endif
1531 tmp1=tmpgrow*phyto(np)-tmpr
1532 tmp2=tmpgrow*phyto(np)*(exp(-sig1*nh4local)+NH4limit(np))
1533 & -tmpr
1534 if (tmp1.ne.0. _d 0) then
1535 Rstarlocal(np)=ksatPO4(np)*tmpr/tmp1
1536 else
1537 Rstarlocal(np)=-9999. _d 0
1538 endif
1539 if (tmp2.ne.0. _d 0) then
1540 RNstarlocal(np)=ksatNO3(np)*
1541 & (tmpr-tmpgrow*NH4limit(np)*phyto(np))/tmp2
1542 else
1543 RNstarlocal(np)=-9999. _d 0
1544 endif
1545 #endif
1546 #ifdef DAR_DIAG_GROW
1547 c include temp, light, nutrients
1548 c Growlocal(np)=grow(np)
1549 c include temp and light, but not nutrients
1550 Growlocal(np)=ngrow(np)*mu(np)*ilimit(np)*
1551 & phytoTempFunction(np)
1552 c include temp, but not nutrients or light
1553 c Growlocal(np)=ngrow(np)*mu(np)*
1554 c & phytoTempFunction(np)
1555 Growsqlocal(np)=Growlocal(np)**2
1556 #endif
1557 enddo
1558 c end np loop
1559 if (debug.eq.10) print*,'dphytodt',dphytodt
1560 c
1561 #ifdef OLD_GRAZE
1562 c zooplankton growth by grazing
1563 do nz=1,nzmax
1564 c zoo in P currency
1565 dzooPdt(nz) = grazingP(nz)*zooP(nz)
1566 C zooplankton stoichiometry varies according to food source
1567 dzooNdt(nz) = grazingN(nz)*zooP(nz)
1568 dzooFedt(nz) = grazingFe(nz)*zooP(nz)
1569 dzooSidt(nz) = grazingSi(nz)*zooP(nz)
1570 enddo
1571 #else
1572 do nz=1,nzmax
1573 c zoo in P currency
1574 dzooPdt(nz) = sumgrazzoo(nz)
1575 C zooplankton stoichiometry varies according to food source
1576 dzooNdt(nz) = sumgrazzooN(nz)
1577 dzooFedt(nz) = sumgrazzooFe(nz)
1578 dzooSidt(nz) = sumgrazzooSi(nz)
1579 enddo
1580 #endif
1581 if (debug.eq.10) print*,'dZooPdt',dZooPdt
1582
1583 c zooplankton mortality
1584 do nz=1,nzmax
1585 c zoo in P currency
1586 dzooPdt(nz) = dzooPdt(nz)
1587 & - mortzoo(nz)*
1588 & mortZTempFunction*zooP(nz)
1589 & - mortzoo2(nz)*
1590 & mortZ2TempFunction*zooP(nz)**2
1591 c zooplankton in other currencies
1592 C zooplankton stoichiometry varies according to food source
1593 dzooNdt(nz) = dzooNdt(nz)
1594 & - mortzoo(nz)*
1595 & mortZTempFunction*zooN(nz)
1596 & - mortzoo2(nz)*
1597 & mortZ2TempFunction*zooN(nz)**2
1598 dzooFedt(nz) = dzooFedt(nz)
1599 & - mortzoo(nz)*
1600 & mortZTempFunction*zooFe(nz)
1601 & - mortzoo2(nz)*
1602 & mortZ2TempFunction*zooFe(nz)**2
1603 dzooSidt(nz) = dzooSidt(nz)
1604 & - mortzoo(nz)*
1605 & mortZTempFunction*zooSi(nz)
1606 & - mortzoo2(nz)*
1607 & mortZ2TempFunction*zooSi(nz)**2
1608 enddo
1609
1610
1611 c sum contributions to inorganic nutrient tendencies
1612 dPO4dt = - consumpPO4 +
1613 #ifdef ALLOW_CDOM
1614 & DOPremin
1615 #else
1616 & preminP + DOPremin
1617 #endif
1618 dNH4dt = - consumpNH4 +
1619 #ifdef ALLOW_CDOM
1620 & DONremin
1621 #else
1622 & preminN + DONremin
1623 #endif
1624 & - NO2prod
1625 dNO2dt = - consumpNO2
1626 & + NO2prod - NO3prod
1627 dNO3dt = - consumpNO3
1628 & + NO3prod
1629 c-ONLYNO3 dNO3dt = - consumpNO3 + preminN + DONremin
1630 #ifdef ALLOW_DENIT
1631 if (O2local.le.O2crit) then
1632 denit = denit_np
1633 #ifdef ALLOW_CDOM
1634 & *(DOPremin)
1635 #else
1636 & *(preminP + DOPremin)
1637 #endif
1638 dNO3dt = dNO3dt - denit
1639 dNH4dt = dNH4dt -
1640 #ifdef ALLOW_CDOM
1641 & (DONremin)
1642 #else
1643 & (preminN + DONremin)
1644 #endif
1645 endif
1646 #endif
1647 dFeTdt = - consumpFeT +
1648 #ifdef ALLOW_CDOM
1649 & DOFeremin
1650 #else
1651 & preminFe + DOFeremin
1652 #endif
1653 #ifdef PART_SCAV
1654 & - scav_part*freefelocal +
1655 #else
1656 & - scav*freefelocal +
1657 #endif
1658 & alpfe*inputFelocal/dzlocal
1659 dSidt = - consumpSi + preminSi
1660
1661 c tendency of dissolved organic pool
1662 dDOPdt = totphy_dop + totzoo_dop - DOPremin
1663 #ifdef ALLOW_CDOM
1664 & +preminP + cdomp_degrd
1665 #endif
1666 dDONdt = totphy_don + totzoo_don - DONremin
1667 #ifdef ALLOW_CDOM
1668 & +preminN + cdomn_degrd
1669 #endif
1670 dDOFedt = totphy_dofe + totzoo_dofe - DOFeremin
1671 #ifdef ALLOW_CDOM
1672 & +preminFe + cdomfe_degrd
1673 #endif
1674
1675 c tendency of particulate detritus pools
1676 dpopdt = totphy_pop + totzoo_pop - preminP + psinkP
1677 #ifdef ALLOW_CDOM
1678 & -preminP_cdom
1679 #endif
1680 dpondt = totphy_pon + totzoo_pon - preminN + psinkN
1681 #ifdef ALLOW_CDOM
1682 & -preminN_cdom
1683 #endif
1684 dpofedt = totphy_pofe + totzoo_pofe - preminFe + psinkFe
1685 #ifdef ALLOW_CDOM
1686 & -preminFe_cdom
1687 #endif
1688 dpSidt = totphy_posi + totzoo_posi - preminSi + psinkSi
1689 #ifdef ALLOW_CARBON
1690 dDICdt = - consumpDIC - consumpDIC_PIC
1691 #ifdef ALLOW_CDOM
1692 & + DOCremin
1693 #else
1694 & + preminC + DOCremin
1695 #endif
1696 & + disscPIC
1697 dDOCdt = totphy_doc + totzoo_doc - DOCremin
1698 #ifdef ALLOW_CDOM
1699 & +preminC + cdomc_degrd
1700 #endif
1701 dPOCdt = totphy_poc + totzoo_poc - preminC + psinkC
1702 #ifdef ALLOW_CDOM
1703 & -preminC_cdom
1704 #endif
1705 dPICdt = totphy_pic + totzoo_pic - disscPIC + psinkPIC
1706 dALKdt = - dNO3dt - 2.d0 * (consumpDIC_PIC - disscPIC)
1707 c should be = O2prod - preminP - DOPremin?
1708 c OLD WAY
1709 c dO2dt = - R_OP*dPO4dt
1710 c production of O2 by photosynthesis
1711 dO2dt = R_OP*consumpPO4
1712 c loss of O2 by remineralization
1713 if (O2local.gt.O2crit) then
1714 dO2dt = dO2dt - R_OP
1715 #ifdef ALLOW_CDOM
1716 & *(DOPremin)
1717 #else
1718 & *(preminP + DOPremin)
1719 #endif
1720 endif
1721 #ifdef OLD_GRAZE
1722 do nz=1,nzmax
1723 dzooCdt(nz) = grazingC(nz)*zooClocal(nz)
1724 & - mortzoo(nz)*
1725 & mortZTempFunction*zooClocal(nz)
1726 & - mortzoo2(nz)*
1727 & mortZ2TempFunction*zooClocal(nz)**2
1728 enddo
1729 #else
1730 do nz=1,nzmax
1731 dzooCdt(nz) = sumgrazzooc(nz)
1732 & - mortzoo(nz)*
1733 & mortZTempFunction*zooClocal(nz)
1734 & - mortzoo2(nz)*
1735 & mortZ2TempFunction*zooClocal(nz)**2
1736 enddo
1737 #endif
1738
1739 #ifdef ALLOW_CDOM
1740 dcdomdt = totphy_cdomp + totzoo_cdomp + preminP_cdom
1741 & -cdomp_degrd
1742 #endif
1743
1744 #endif
1745
1746 if (debug.eq.10) print*,'dDOPdt', dDOPdt
1747 if (debug.eq.10) print*,'dpopdt',dpopdt
1748 if (debug.eq.10) print*,'dDONdt',dDONdt
1749 if (debug.eq.10) print*,'dpondt',dpondt
1750 c
1751 c -------------------------------------------------------------------
1752 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1753 c --------------------------------------------------------------------------
1754
1755 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-
1756 c Mutation - apply mutation to tendencies [jbmodif]
1757
1758 #ifdef ALLOW_MUTANTS
1759 c apply to all sisters when first sister is encountered
1760 if(runtim .gt. threeyr) then
1761 mutfor=1 _d -8
1762 mutback=1 _d -12
1763 if(numtax .gt. 1)then
1764 do np=1,npro
1765 if(mod(np,numtax).eq. 1. _d 0)then
1766 nsisone = np
1767 nsistwo = np+1
1768 nsisthree = np+2
1769 nsisfour = np+3
1770
1771 grow1 = PspecificPO4(nsisone)
1772 grow2 = PspecificPO4(nsistwo)
1773
1774 if(numtax.eq.2)grow3 = 0.0 _d 0
1775 if(numtax.eq.2)grow4 = 0.0 _d 0
1776
1777 if(numtax.eq.3)grow4 = 0.0 _d 0
1778 if(numtax.ge.3)grow3 = PspecificPO4(nsisthree)
1779
1780 if(numtax.eq.4)grow4 = PspecificPO4(nsisfour)
1781
1782
1783
1784 dphytodt(nsisone) = dphytodt(nsisone)
1785 & - grow1 *1.4427 _d 0*mutfor
1786 & - grow1 *1.4427 _d 0*mutfor
1787 & - grow1 *1.4427 _d 0*mutfor
1788 & + grow2 *1.4427 _d 0*mutback
1789 & + grow3 *1.4427 _d 0*mutback
1790 & + grow4 *1.4427 _d 0*mutback
1791
1792 dphytodt(nsistwo) = dphytodt(nsistwo)
1793 & - grow2 *1.4427 _d 0*mutback
1794 & + grow1 *1.4427 _d 0*mutfor
1795
1796 if(numtax .ge. 3)then
1797 dphytodt(nsisthree) = dphytodt(nsisthree)
1798 & - grow3 *1.4427 _d 0*mutback
1799 & + grow1 *1.4427 _d 0*mutfor
1800 endif
1801
1802 if(numtax .eq. 4)then
1803 dphytodt(nsisfour) = dphytodt(nsisfour)
1804 & - grow4 *1.4427 _d 0*mutback
1805 & + grow1 *1.4427 _d 0*mutfor
1806 c QQQQQQQQQQ FIX FOR NIT RUNS ONLY!!!
1807 if (phyto(nsisfour).eq.0. _d 0) then
1808 if (phyto(nsistwo).eq.0. _d 0) then
1809 if (dphytodt(nsistwo).gt.dphytodt(nsisfour)) then
1810 dphytodt(nsisfour)=dphytodt(nsistwo)
1811 endif
1812 endif
1813 if (phyto(nsisthree).eq.0. _d 0) then
1814 if (dphytodt(nsisthree).gt.dphytodt(nsisfour)) then
1815 dphytodt(nsisfour)=dphytodt(nsisthree)
1816 endif
1817 endif
1818 endif
1819 c QQQQQQQQQQQQQ
1820 endif
1821
1822 c QQQQQQQQQQQQTEST
1823 if (debug.eq.11) then
1824 if (PARlocal.gt.1. _d 0) then
1825 if (dphytodt(nsistwo).gt.dphytodt(nsisfour).and.
1826 & dphytodt(nsisfour).gt.0. _d 0) then
1827 print*,'QQQQ nsistwo>nsisfour',nsistwo,nsisfour,
1828 & dphytodt(nsistwo), dphytodt(nsisfour),
1829 & phyto(nsistwo), phyto(nsisfour),
1830 & phyto(nsisone)
1831 endif
1832 if (dphytodt(nsisthree).gt.dphytodt(nsisfour).and.
1833 & dphytodt(nsisfour).gt.0. _d 0) then
1834 print*,'QQQQ nsisthree>nsisfour',nsisthree,nsisfour,
1835 & dphytodt(nsisthree), dphytodt(nsisfour),
1836 & phyto(nsisthree), phyto(nsisfour),
1837 & phyto(nsisone)
1838 endif
1839 if (dphytodt(nsisfour).gt.dphytodt(nsisone).and.
1840 & dphytodt(nsisone).gt.0. _d 0) then
1841 print*,' BIG QQQQ nsisfour>nsisone',nsisone,nsisfour,
1842 & dphytodt(nsisfour), dphytodt(nsisone),
1843 & phyto(nsisfour), phyto(nsisone)
1844 endif
1845 endif
1846 endif
1847 c QQQQQQQQQTEST
1848 endif
1849 enddo
1850 endif
1851 endif
1852
1853 c mutation is finished
1854 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-
1855 #endif
1856
1857
1858
1859 RETURN
1860 END
1861 #endif /*MONOD*/
1862 #endif /*ALLOW_PTRACERS*/
1863 c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22