/[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.4 - (show annotations) (download)
Fri Oct 21 20:56:53 2011 UTC (13 years, 9 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63d_20111107
Changes since 1.3: +30 -12 lines
o add quadratic mortality for zooplankton

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

  ViewVC Help
Powered by ViewVC 1.1.22