/[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.1 - (show annotations) (download)
Wed Apr 13 18:56:25 2011 UTC (14 years, 4 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_baseline, ctrb_darwin2_ckpt62z_20110622
darwin2 initial checkin

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

  ViewVC Help
Powered by ViewVC 1.1.22