/[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.5 - (show annotations) (download)
Tue Nov 15 20:16:34 2011 UTC (13 years, 8 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63k_20120317
Changes since 1.4: +1 -2 lines
remove unnused header

1 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F,v 1.4 2011/10/21 20:56:53 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 "MONOD_SIZE.h"
98 #include "MONOD.h"
99 #include "DARWIN_PARAMS.h"
100
101 c ANNA set wavebands params
102 #ifdef WAVEBANDS
103 #include "SPECTRAL_SIZE.h"
104 #include "WAVEBANDS_PARAMS.h"
105 #endif
106
107
108 C !INPUT PARAMETERS: ===================================================
109 C myThid :: thread number
110 INTEGER myThid
111 CEOP
112 c === GLOBAL VARIABLES =====================
113 c npmax = no of phyto functional groups
114 c nzmax = no of grazer species
115 c phyto = phytoplankton
116 c zoo = zooplankton
117 _RL phyto(npmax)
118 _RL zooP(nzmax)
119 _RL zooN(nzmax)
120 _RL zooFe(nzmax)
121 _RL zooSi(nzmax)
122 _RL PP
123 _RL Nfix
124 _RL denit
125 _RL Chl
126 _RL PO4local
127 _RL NO3local
128 _RL FeTlocal
129 _RL Silocal
130 _RL NO2local
131 _RL NH4local
132 _RL DOPlocal
133 _RL DONlocal
134 _RL DOFelocal
135 _RL POPlocal
136 _RL PONlocal
137 _RL POFelocal
138 _RL PSilocal
139 _RL phytoup(npmax)
140 _RL POPuplocal
141 _RL PONuplocal
142 _RL POFeuplocal
143 _RL PSiuplocal
144 _RL PARlocal
145 _RL Tlocal
146 _RL Slocal
147 _RL freefelocal
148 _RL inputFelocal
149 _RL bottom
150 _RL dzlocal
151 _RL Rstarlocal(npmax)
152 _RL RNstarlocal(npmax)
153 #ifdef DAR_DIAG_GROW
154 _RL Growlocal(npmax)
155 _RL Growsqlocal(npmax)
156 #endif
157 #ifdef ALLOW_DIAZ
158 #ifdef DAR_DIAG_NFIXP
159 _RL NfixPlocal(npmax)
160 #endif
161 #endif
162 INTEGER debug
163 _RL dphytodt(npmax)
164 _RL dzooPdt(nzmax)
165 _RL dzooNdt(nzmax)
166 _RL dzooFedt(nzmax)
167 _RL dzooSidt(nzmax)
168 _RL dPO4dt
169 _RL dNO3dt
170 _RL dNO2dt
171 _RL dNH4dt
172 _RL dFeTdt
173 _RL dSidt
174 _RL dDOPdt
175 _RL dDONdt
176 _RL dDOFedt
177 _RL dPOPdt
178 _RL dPONdt
179 _RL dPOFedt
180 _RL dPSidt
181 #ifdef ALLOW_CARBON
182 _RL DIClocal
183 _RL DOClocal
184 _RL POClocal
185 _RL PIClocal
186 _RL ALKlocal
187 _RL O2local
188 _RL ZooClocal(nzmax)
189 _RL POCuplocal
190 _RL PICuplocal
191 _RL dDICdt
192 _RL dDOCdt
193 _RL dPOCdt
194 _RL dPICdt
195 _RL dALKdt
196 _RL dO2dt
197 _RL dZOOCdt(nzmax)
198 #endif
199 #ifdef GEIDER
200 _RL phychl(npmax)
201 #ifdef DYNAMIC_CHL
202 _RL dphychl(npmax)
203 _RL Chlup(npmax)
204 #endif
205 #endif
206 #ifdef ALLOW_PAR_DAY
207 _RL PARdaylocal
208 #endif
209 #ifdef DAR_DIAG_CHL
210 _RL ChlGeiderlocal, ChlDoneylocal, ChlCloernlocal
211 #endif
212 _RL runtim
213
214 c ANNA Global variables for WAVEBANDS
215 c ANNA these variables are passed in/out of darwin_forcing.F
216 #ifdef WAVEBANDS
217 _RL PARwlocal(tlam) !PAR at midpoint of previous(in) and local(out) gridcell
218 #endif
219 c ANNA endif
220
221
222
223
224
225 c LOCAL VARIABLES
226 c -------------------------------------------------------------
227
228 c WORKING VARIABLES
229 c np = phytoplankton index
230 integer np
231 c nz = zooplankton index
232 integer nz
233
234 c variables for phytoplankton growth rate/nutrient limitation
235 c phytoplankton specific nutrient limitation term
236 _RL limit(npmax)
237 c phytoplankton light limitation term
238 _RL ilimit(npmax)
239 _RL ngrow(npmax)
240 _RL grow(npmax)
241 _RL PspecificPO4(npmax)
242 _RL phytoTempFunction(npmax)
243 _RL dummy
244 _RL Ndummy
245 _RL Nsourcelimit(npmax)
246 _RL Nlimit(npmax)
247 _RL NO3limit(npmax)
248 _RL NO2limit(npmax)
249 _RL NH4limit(npmax)
250
251 c for check N limit scheme
252 _RL Ndiff
253 _RL NO3limcheck
254 _RL NO2limcheck
255 _RL Ndummy1
256 LOGICAL check_nlim
257
258 #ifndef OLD_NSCHEME
259 c [jbmodif] some new N terms
260 integer N2only
261 integer noNOdadv
262 integer NOreducost
263 _RL NO2zoNH4
264 _RL NOXzoNH4
265 #endif
266
267 c varible for mimumum phyto
268 _RL phytomin(npmax)
269
270 #ifdef OLD_GRAZE
271 c variables for zooplankton grazing rates
272 _RL zooTempFunction(nzmax)
273 _RL grazing_phyto(npmax)
274 _RL grazingP(nzmax)
275 _RL grazingN(nzmax)
276 _RL grazingFe(nzmax)
277 _RL grazingSi(nzmax)
278 #else
279 c variables for zooplankton grazing rates
280 _RL zooTempFunction(nzmax)
281 _RL allphyto(nzmax)
282 _RL denphyto(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 #ifdef SER_GRAZ
955 denphyto(nz)=0. _d 0
956 do np=1,npmax
957 denphyto(nz)=denphyto(nz)+
958 & (palat(np,nz)*phyto(np)/allphyto(nz))*phyto(np)
959 enddo
960 if (denphyto(nz).le.0. _d 0) denphyto(nz)=phygrazmin
961 #endif
962 do np=1,npmax
963 tmpz=max(0. _d 0,(allphyto(nz)-phygrazmin) )
964 grazphy(np,nz)=grazemax(nz)*
965 #ifdef SER_GRAZ
966 c as in Vallina et al, 2011
967 & (((palat(np,nz)*phyto(np)/allphyto(nz))*phyto(np))/
968 & denphyto(nz)) *
969 #else
970 c as in Dutkiewicz et al, GBC, 2009
971 & (palat(np,nz)*phyto(np)/allphyto(nz))*
972 #endif
973 & ( tmpz**hollexp/
974 & (tmpz**hollexp+kgrazesat**hollexp) )
975 enddo
976 enddo
977 if (debug.eq.2) print*,'allphyto',allphyto
978 c if (debug.eq.2) print*,'grazephy',grazphy
979 c sum over zoo for impact on phyto
980 do np=1,npmax
981 sumgrazphy(np)=0. _d 0
982 do nz=1,nzmax
983 sumgrazphy(np)=sumgrazphy(np)+
984 & grazphy(np,nz)*zooP(nz)
985 enddo
986 enddo
987 if (debug.eq.2) print*,'sumgrazephy',sumgrazphy
988 c sum over phy for impact on zoo, and all remainder to go to POM
989 do nz=1,nzmax
990 sumgrazzoo(nz)=0. _d 0
991 sumgrazzooN(nz)=0. _d 0
992 sumgrazzooFe(nz)=0. _d 0
993 sumgrazzooSi(nz)=0. _d 0
994 sumgrazloss(nz)=0. _d 0
995 sumgrazlossN(nz)=0. _d 0
996 sumgrazlossFe(nz)=0. _d 0
997 sumgrazlossSi(nz)=0. _d 0
998 #ifdef ALLOW_CARBON
999 sumgrazzooC(nz)=0. _d 0
1000 sumgrazlossC(nz)=0. _d 0
1001 sumgrazlossPIC(nz)=0. _d 0
1002 #endif
1003 do np=1,npmax
1004 sumgrazzoo(nz)=sumgrazzoo(nz)+
1005 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)
1006 sumgrazloss(nz)=sumgrazloss(nz)+
1007 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*zooP(nz)
1008 sumgrazzooN(nz)=sumgrazzooN(nz)+
1009 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_NP(np)
1010 sumgrazlossN(nz)=sumgrazlossN(nz)+
1011 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1012 & zooP(nz)*R_NP(np)
1013 sumgrazzooFe(nz)=sumgrazzooFe(nz)+
1014 & asseff(np,nz)*grazphy(np,nz)*
1015 & zooP(nz)*R_FeP(np)
1016 sumgrazlossFe(nz)=sumgrazlossFe(nz)+
1017 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1018 & zooP(nz)*R_FeP(np)
1019 sumgrazzooSi(nz)=sumgrazzooSi(nz)+
1020 & asseff(np,nz)*grazphy(np,nz)*
1021 & zooP(nz)*R_SiP(np)
1022 sumgrazlossSi(nz)=sumgrazlossSi(nz)+
1023 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1024 & zooP(nz)*R_SiP(np)
1025 #ifdef ALLOW_CARBON
1026 sumgrazzooC(nz)=sumgrazzooC(nz)+
1027 & asseff(np,nz)*grazphy(np,nz)*zooP(nz)*R_PC(np)
1028 sumgrazlossC(nz)=sumgrazlossC(nz)+
1029 & (1. _d 0-asseff(np,nz))*grazphy(np,nz)*
1030 & zooP(nz)*R_PC(np)
1031 sumgrazlossPIC(nz)=sumgrazlossPIC(nz)+
1032 & (1. _d 0)*grazphy(np,nz)*
1033 & zooP(nz)*R_PC(np)*R_PICPOC(np)
1034 #endif
1035 enddo
1036 enddo
1037 if (debug.eq.2) print*,'sumgrazzoo',sumgrazzoo
1038 if (debug.eq.2) print*,'sumgrazloss',sumgrazloss
1039 if (debug.eq.2) print*,'sumgrazzooN',sumgrazzooN
1040 if (debug.eq.2) print*,'sumgrazlossN',sumgrazlossN
1041 if (debug.eq.2) print*,'sumgrazzooFe',sumgrazzooFe
1042 if (debug.eq.2) print*,'sumgrazlossFe',sumgrazlossFe
1043 if (debug.eq.2) print*,'sumgrazzooSi',sumgrazzooSi
1044 if (debug.eq.2) print*,'sumgrazlossSi',sumgrazlossSi
1045 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1046 #endif
1047
1048 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1049 c accumulate particulate and dissolved detritus
1050 do np=1, npmax
1051 totphy_pop=totphy_pop+
1052 & ExportFracP(np)*mortphy(np)*phytomin(np)
1053 totphy_dop=totphy_dop+
1054 & (1. _d 0-ExportFracP(np))*mortphy(np)*phytomin(np)
1055 totphy_pon=totphy_pon+ R_NP(np)*
1056 & ExportFracP(np)*mortphy(np)*phytomin(np)
1057 totphy_don=totphy_don+ R_NP(np)*
1058 & (1. _d 0-ExportFracP(np))*mortphy(np)*phytomin(np)
1059 totphy_pofe=totphy_pofe+ R_FeP(np)*
1060 & ExportFracP(np)*mortphy(np)*phytomin(np)
1061 totphy_dofe=totphy_dofe+ R_FeP(np)*
1062 & (1. _d 0-ExportFracP(np))*mortphy(np)*phytomin(np)
1063 totphy_posi=totphy_posi+ R_SiP(np)*
1064 & mortphy(np)*phytomin(np)
1065 #ifdef ALLOW_CARBON
1066 totphy_poc=totphy_poc+ R_PC(np)*
1067 & ExportFracP(np)*mortphy(np)*phytomin(np)
1068 totphy_doc=totphy_doc+ R_PC(np)*
1069 & (1. _d 0-ExportFracP(np))*mortphy(np)*phytomin(np)
1070 totphy_pic=totphy_pic+ R_PC(np)*R_PICPOC(np)*
1071 & mortphy(np)*phytomin(np)
1072 #endif
1073 enddo
1074 if (debug.eq.3) print*,'tot_phy_pop',totphy_pop
1075 if (debug.eq.3) print*,'tot_phy_dop',totphy_dop
1076 if (debug.eq.3) print*,'tot_phy_pon',totphy_pon
1077 if (debug.eq.3) print*,'tot_phy_don',totphy_don
1078 if (debug.eq.3) print*,'tot_phy_pofe',totphy_pofe
1079 if (debug.eq.3) print*,'tot_phy_dofe',totphy_dofe
1080 if (debug.eq.3) print*,'tot_phy_posi',totphy_posi
1081
1082
1083 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1084
1085
1086 #ifdef OLD_GRAZE
1087 c ****************** ZOO GRAZING RATE ****************************
1088 c determine zooplankton grazing rates
1089 do nz = 1, nzmax
1090 c grazing: sum contribution from all phytoplankton
1091 grazingP(nz) = 0.0 _d 0
1092 grazingN(nz) = 0.0 _d 0
1093 grazingFe(nz) = 0.0 _d 0
1094 grazingSi(nz) = 0.0 _d 0
1095 #ifdef ALLOW_CARBON
1096 grazingC(nz) = 0.0 _d 0
1097 #endif
1098 do np = 1, npmax
1099 facpz = (phytomin(np)/(phytomin(np) + kgrazesat))
1100 & *zooTempFunction(nz)
1101 grazingP(nz) = grazingP(nz) +
1102 & graze(np,nz)*facpz
1103 grazingN(nz) = grazingN(nz) +
1104 & graze(np,nz)*R_NP(np)*facpz
1105 grazingFe(nz) = grazingFe(nz) +
1106 & graze(np,nz)*R_FeP(np)*facpz
1107 grazingSi(nz) = grazingSi(nz) +
1108 & graze(np,nz)*R_SiP(np)*facpz
1109 #ifdef ALLOW_CARBON
1110 grazingC(nz) = grazingC(nz) +
1111 & graze(np,nz)*R_PC(np)*facpz
1112 #endif
1113 enddo
1114 enddo
1115 if (debug.eq.4) print*,'grazingP', grazingP
1116 if (debug.eq.4) print*,'grazingN', grazingN
1117 if (debug.eq.4) print*,'grazingFe', grazingFe
1118 if (debug.eq.4) print*,'grazingSi', grazingSi
1119 c ************* END ZOO GRAZING *********************************
1120 #endif
1121 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1122 c accumulate particulate and dissolved detritus
1123 do nz=1, nzmax
1124 totzoo_pop=totzoo_pop+
1125 & ExportFracZ(nz)*( mortzoo(nz)*zooP(nz)
1126 & + mortzoo2(nz)*zooP(nz)**2 )
1127 totzoo_dop=totzoo_dop+
1128 & (1. _d 0-ExportFracZ(nz))*(
1129 & mortzoo(nz)*zooP(nz)+
1130 & mortzoo2(nz)*zooP(nz)**2 )
1131 totzoo_pon=totzoo_pon+
1132 & ExportFracZ(nz)*( mortzoo(nz)*zooN(nz)
1133 & + mortzoo2(nz)*zooN(nz)**2 )
1134 totzoo_don=totzoo_don+
1135 & (1. _d 0-ExportFracZ(nz))*(
1136 & mortzoo(nz)*zooN(nz)+
1137 & mortzoo2(nz)*zooN(nz)**2 )
1138 totzoo_pofe=totzoo_pofe+
1139 & ExportFracZ(nz)*( mortzoo(nz)*zooFe(nz)
1140 & + mortzoo2(nz)*zooFe(nz)**2 )
1141 totzoo_dofe=totzoo_dofe+
1142 & (1. _d 0-ExportFracZ(nz))*(
1143 & mortzoo(nz)*zooFe(nz) +
1144 & mortzoo2(nz)*zooFe(nz)**2 )
1145 totzoo_posi=totzoo_posi+
1146 & ( mortzoo(nz)*zooSi(nz)+
1147 & mortzoo2(nz)*zooSi(nz)**2 )
1148 #ifdef ALLOW_CARBON
1149 totzoo_poc=totzoo_poc+
1150 & ExportFracZ(nz)*( mortzoo(nz)*zooClocal(nz)
1151 & + mortzoo2(nz)*zooClocal(nz)**2 )
1152 totzoo_doc=totzoo_doc+
1153 & (1. _d 0-ExportFracZ(nz))*( mortzoo(nz)*zooClocal(nz)
1154 & + mortzoo2(nz)*zooClocal(nz)**2 )
1155 #endif
1156 enddo
1157
1158 #ifndef OLD_GRAZE
1159 do nz=1, nzmax
1160 totzoo_pop=totzoo_pop+
1161 & ExportFracGraz(nz)*sumgrazloss(nz)
1162 totzoo_dop=totzoo_dop+
1163 & (1. _d 0-ExportFracGraz(nz))*sumgrazloss(nz)
1164 totzoo_pon=totzoo_pon+
1165 & ExportFracGraz(nz)*sumgrazlossN(nz)
1166 totzoo_don=totzoo_don+
1167 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossN(nz)
1168 totzoo_pofe=totzoo_pofe+
1169 & ExportFracGraz(nz)*sumgrazlossFe(nz)
1170 totzoo_dofe=totzoo_dofe+
1171 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossFe(nz)
1172 totzoo_posi=totzoo_posi+
1173 & sumgrazlossSi(nz)
1174 #ifdef ALLOW_CARBON
1175 totzoo_poc=totzoo_poc+
1176 & ExportFracGraz(nz)*sumgrazlossC(nz)
1177 totzoo_doc=totzoo_doc+
1178 & (1. _d 0-ExportFracGraz(nz))*sumgrazlossC(nz)
1179 totzoo_pic=totzoo_pic+
1180 & sumgrazlossPIC(nz)
1181 #endif
1182 enddo
1183 #endif
1184 if (debug.eq.5) print*,'totzoo_pop',totzoo_pop
1185 if (debug.eq.5) print*,'totzoo_dop',totzoo_dop
1186 if (debug.eq.5) print*,'totzoo_pon',totzoo_pon
1187 if (debug.eq.5) print*,'totzoo_don',totzoo_don
1188 if (debug.eq.5) print*,'totzoo_pofe',totzoo_pofe
1189 if (debug.eq.5) print*,'totzoo_dofe',totzoo_dofe
1190 if (debug.eq.5) print*,'totzoo_posi',totzoo_posi
1191 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1192
1193 c ********************* NUTRIENT UPTAKE *******************************
1194 c determine nutrient uptake
1195 c consumption - sum of phytoplankton contributions
1196 do np = 1, npmax
1197 c phospate uptake by each phytoplankton
1198 #ifndef GEIDER
1199 grow(np)=ngrow(np)*mu(np)*limit(np)*ilimit(np)*
1200 & phytoTempFunction(np)
1201 #endif
1202 #ifdef GEIDER
1203 grow(np)=ngrow(np)*pcarbon(np)
1204 if (debug.eq.1) print*,'grow', grow(np), pcarbon(np)
1205 if (debug.eq.14) print*,'grow', grow(np), pcarbon(np)
1206 #ifdef DYNAMIC_CHL
1207 c geider 97 for dChl/dt (source part) Eq. 3
1208 if (acclim(np).gt. 0. _d 0) then
1209 rhochl(np)=chl2cmax(np) *
1210 & (grow(np)/(alpha_I(np)*acclim(np)) )
1211 else
1212 rhochl(np)= 0. _d 0
1213 endif
1214 if (debug.eq.14) print*,'rhochl',rhochl(np)
1215 #endif
1216 #endif
1217 PspecificPO4(np) = grow(np)*phyto(np)
1218 c write(6,*)'np =',np, ' PspecificPO4 ='
1219 c & ,PspecificPO4(np)
1220 consumpPO4 = consumpPO4 + PspecificPO4(np)
1221 consumpFeT = consumpFeT + PspecificPO4(np)*R_FeP(np)
1222 consumpSi = consumpSi + PspecificPO4(np)*R_SiP(np)
1223 cswd should have O2prod as function of np?
1224 c New Way of doing Nitrogen Consumption .......................
1225 if(diazotroph(np) .ne. 1.0 _d 0)then
1226 if (Nlimit(np).le.0. _d 0) then
1227 consumpNO3 = consumpNO3
1228 consumpNO2 = consumpNO2
1229 consumpNH4 = consumpNH4
1230 else
1231 consumpNO3 = consumpNO3 +
1232 & NO3limit(np)/Nlimit(np)*PspecificPO4(np)*R_NP(np)
1233 consumpNO2 = consumpNO2 +
1234 & NO2limit(np)/Nlimit(np)* PspecificPO4(np)*R_NP(np)
1235 consumpNH4 = consumpNH4 +
1236 & NH4limit(np)/Nlimit(np)*PspecificPO4(np)*R_NP(np)
1237 endif
1238 else
1239 consumpNO3 = consumpNO3
1240 consumpNO2 = consumpNO2
1241 consumpNH4 = consumpNH4
1242 Nfix=Nfix+PspecificPO4(np)*R_NP(np)
1243 #ifdef ALLOW_DIAZ
1244 #ifdef DAR_DIAG_NFIXP
1245 NfixPlocal(np)=PspecificPO4(np)*R_NP(np)
1246 #endif
1247 #endif
1248 endif
1249 #ifdef ALLOW_CARBON
1250 consumpDIC = consumpDIC + PspecificPO4(np)*R_PC(np)
1251 consumpDIC_PIC = consumpDIC_PIC +
1252 & PspecificPO4(np)*R_PC(np)*R_PICPOC(np)
1253 #endif
1254 enddo
1255 if (debug.eq.7) print*,'local', parlocal,tlocal,po4local,
1256 & no3local, no2local,nh4local,fetlocal,silocal
1257 if (debug.eq.7) print*,'grow',grow
1258 if (debug.eq.6) print*,'pspecificpo4', PspecificPO4
1259 if (debug.eq.6) print*,'consumpPO4', consumpPO4
1260 if (debug.eq.6) print*,'consumpFeT', consumpFeT
1261 if (debug.eq.6) print*,'consumpSi ', consumpsi
1262 if (debug.eq.6) print*,'consumpNO3', consumpNO3
1263 if (debug.eq.6) print*,'consumpNO2', consumpNO2
1264 if (debug.eq.6) print*,'consumpNH4', consumpNH4
1265 c ****************** END NUTRIENT UPTAKE ****************************
1266
1267 c sinking phytoplankton and POM
1268 if(bottom .eq. 1.0 _d 0)then
1269 psinkP = (wp_sink*POPuplocal)/(dzlocal)
1270 psinkN = (wn_sink*PONuplocal)/(dzlocal)
1271 psinkFe = (wfe_sink*POFeuplocal)/(dzlocal)
1272 psinkSi = (wsi_sink*PSiuplocal)/(dzlocal)
1273 do np=1,npmax
1274 psinkPhy(np) =
1275 & (wsink(np)*Phytoup(np))/(dzlocal)
1276 enddo
1277 #ifdef DYNAMIC_CHL
1278 do np=1,npmax
1279 psinkChl(np) =
1280 & (wsink(np)*Chlup(np))/(dzlocal)
1281 enddo
1282 #endif
1283 #ifdef ALLOW_CARBON
1284 psinkC = (wc_sink*POCuplocal)/(dzlocal)
1285 psinkPIC = (wpic_sink*PICuplocal)/(dzlocal)
1286 #endif
1287 else
1288 psinkP = (wp_sink*(POPuplocal-POPlocal))/(dzlocal)
1289 psinkN = (wn_sink*(PONuplocal-PONlocal))/(dzlocal)
1290 psinkFe = (wfe_sink*(POFeuplocal-POFelocal))/(dzlocal)
1291 psinkSi = (wsi_sink*(PSiuplocal-PSilocal))/(dzlocal)
1292 do np=1,npmax
1293 psinkPhy(np) =
1294 & (wsink(np))*(Phytoup(np)-Phyto(np))/(dzlocal)
1295 enddo
1296 #ifdef DYNAMIC_CHL
1297 do np=1,npmax
1298 psinkChl(np) =
1299 & (wsink(np))*(Chlup(np)-phychl(np))/(dzlocal)
1300 enddo
1301 #endif
1302 #ifdef ALLOW_CARBON
1303 psinkC = (wc_sink*(POCuplocal-POClocal))/(dzlocal)
1304 psinkPIC = (wpic_sink*(PICuplocal-PIClocal))/(dzlocal)
1305 #endif
1306 endif
1307
1308 c DOM remineralization rates
1309 DOPremin = reminTempFunction * Kdop * DOPlocal
1310 DONremin = reminTempFunction * Kdon * DONlocal
1311 DOFeremin = reminTempFunction * KdoFe * DOFelocal
1312
1313 c remineralization of sinking particulate
1314 preminP = reminTempFunction * Kpremin_P*POPlocal
1315 preminN = reminTempFunction * Kpremin_N*PONlocal
1316 preminFe = reminTempFunction * Kpremin_Fe*POFelocal
1317 preminSi = reminTempFunction * Kpremin_Si*PSilocal
1318
1319 #ifdef ALLOW_CARBON
1320 DOCremin = reminTempFunction * Kdoc * DOClocal
1321 preminC = reminTempFunction * Kpremin_C*POClocal
1322 c dissolution
1323 disscPIC = Kdissc*PIClocal
1324 #endif
1325
1326 c chemistry
1327 c NH4 -> NO2 -> NO3 by bacterial action
1328 NO2prod = knita*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1329 & *NH4local
1330 NO3prod = knitb*( 1. _d 0-min(PARlocal/PAR0,1. _d 0) )
1331 & *NO2local
1332 c NO2prod = knita*NH4local
1333 c NO3prod = knitb*NO2local
1334 c
1335 #ifdef PART_SCAV
1336 scav_poc=POPlocal/1.1321 _d -4
1337 c scav rate
1338 scav_part=scav_rat*scav_inter*(scav_poc**scav_exp)
1339 #endif
1340 c -------------------------------------------------------------------
1341 c calculate tendency terms (and some diagnostics)
1342 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1343 c phytoplankton
1344 do np=1,npmax
1345 dphytodt(np) = PspecificPO4(np)
1346 #ifdef OLD_GRAZE
1347 & - grazing_phyto(np)*
1348 & (phytomin(np)/(phytomin(np) + kgrazesat))
1349 #else
1350 & - sumgrazphy(np)
1351 #endif
1352 & - mortphy(np)*phytomin(np)
1353 & + psinkphy(np)
1354 #ifdef GEIDER
1355 #ifdef DYNAMIC_CHL
1356 dphychl(np) = acclim(np)*PspecificPO4(np)*R_PC(np)
1357 c dphychl(np) = rhochl(np)*PspecificPO4(np)*R_PC(np)
1358 & + acclimtimescl *
1359 & (acclim(np)-chl2c(np))*phyto(np)*R_PC(np)
1360 & +(
1361 #ifdef OLD_GRAZE
1362 & - grazing_phyto(np)*
1363 & (phytomin(np)/(phytomin(np) + kgrazesat))
1364 #else
1365 & - sumgrazphy(np)
1366 #endif
1367 & - mortphy(np)*phytomin(np)) *chl2c(np)*R_PC(np)
1368 & + psinkChl(np)
1369 #endif
1370 Chl=Chl + phychl(np)
1371 #endif
1372 c %% diagnostics
1373 PP = PP + PspecificPO4(np)
1374 c%%%
1375 #ifdef OLD_GRAZE
1376 tmpr=grazing_phyto(np)*
1377 & (phytomin(np)/(phytomin(np) + kgrazesat))
1378 & + mortphy(np)*phytomin(np)
1379 & - psinkphy(np)
1380 #else
1381 tmpr=sumgrazphy(np)
1382 & + mortphy(np)*phytomin(np)
1383 & - psinkphy(np)
1384 #endif
1385 #ifdef DAR_DIAG_RSTAR
1386 #ifndef GEIDER
1387 tmpgrow=ngrow(np)*mu(np)*ilimit(np)*
1388 & phytoTempFunction(np)
1389 #endif
1390 #ifdef GEIDER
1391 tmpgrow=grow(np)/limit(np)
1392 #endif
1393 tmp1=tmpgrow*phyto(np)-tmpr
1394 tmp2=tmpgrow*phyto(np)*(exp(-sig1*nh4local)+NH4limit(np))
1395 & -tmpr
1396 if (tmp1.ne.0. _d 0) then
1397 Rstarlocal(np)=ksatPO4(np)*tmpr/tmp1
1398 else
1399 Rstarlocal(np)=-9999. _d 0
1400 endif
1401 if (tmp2.ne.0. _d 0) then
1402 RNstarlocal(np)=ksatNO3(np)*
1403 & (tmpr-tmpgrow*NH4limit(np)*phyto(np))/tmp2
1404 else
1405 RNstarlocal(np)=-9999. _d 0
1406 endif
1407 #endif
1408 #ifdef DAR_DIAG_GROW
1409 c include temp, light, nutrients
1410 c Growlocal(np)=grow(np)
1411 c include temp and light, but not nutrients
1412 Growlocal(np)=ngrow(np)*mu(np)*ilimit(np)*
1413 & phytoTempFunction(np)
1414 c include temp, but not nutrients or light
1415 c Growlocal(np)=ngrow(np)*mu(np)*
1416 c & phytoTempFunction(np)
1417 Growsqlocal(np)=Growlocal(np)**2
1418 #endif
1419 enddo
1420 c end np loop
1421 if (debug.eq.10) print*,'dphytodt',dphytodt
1422 c
1423 #ifdef OLD_GRAZE
1424 c zooplankton growth by grazing
1425 do nz=1,nzmax
1426 c zoo in P currency
1427 dzooPdt(nz) = grazingP(nz)*zooP(nz)
1428 C zooplankton stoichiometry varies according to food source
1429 dzooNdt(nz) = grazingN(nz)*zooP(nz)
1430 dzooFedt(nz) = grazingFe(nz)*zooP(nz)
1431 dzooSidt(nz) = grazingSi(nz)*zooP(nz)
1432 enddo
1433 #else
1434 do nz=1,nzmax
1435 c zoo in P currency
1436 dzooPdt(nz) = sumgrazzoo(nz)
1437 C zooplankton stoichiometry varies according to food source
1438 dzooNdt(nz) = sumgrazzooN(nz)
1439 dzooFedt(nz) = sumgrazzooFe(nz)
1440 dzooSidt(nz) = sumgrazzooSi(nz)
1441 enddo
1442 #endif
1443 if (debug.eq.10) print*,'dZooPdt',dZooPdt
1444
1445 c zooplankton mortality
1446 do nz=1,nzmax
1447 c zoo in P currency
1448 dzooPdt(nz) = dzooPdt(nz)
1449 & - mortzoo(nz)*zooP(nz)
1450 & - mortzoo2(nz)*zooP(nz)**2
1451 c zooplankton in other currencies
1452 C zooplankton stoichiometry varies according to food source
1453 dzooNdt(nz) = dzooNdt(nz)
1454 & - mortzoo(nz)*zooN(nz)
1455 & - mortzoo2(nz)*zooN(nz)**2
1456 dzooFedt(nz) = dzooFedt(nz)
1457 & - mortzoo(nz)*zooFe(nz)
1458 & - mortzoo2(nz)*zooFe(nz)**2
1459 dzooSidt(nz) = dzooSidt(nz)
1460 & - mortzoo(nz)*zooSi(nz)
1461 & - mortzoo2(nz)*zooSi(nz)**2
1462 enddo
1463
1464
1465 c sum contributions to inorganic nutrient tendencies
1466 dPO4dt = - consumpPO4 + preminP + DOPremin
1467 dNH4dt = - consumpNH4 + preminN + DONremin
1468 & - NO2prod
1469 dNO2dt = - consumpNO2
1470 & + NO2prod - NO3prod
1471 dNO3dt = - consumpNO3
1472 & + NO3prod
1473 c-ONLYNO3 dNO3dt = - consumpNO3 + preminN + DONremin
1474 #ifdef ALLOW_DENIT
1475 if (O2local.le.O2crit) then
1476 denit = denit_np*(preminP + DOPremin)
1477 dNO3dt = dNO3dt - denit
1478 dNH4dt = dNH4dt - (preminN + DONremin)
1479 endif
1480 #endif
1481 dFeTdt = - consumpFeT + preminFe + DOFeremin
1482 #ifdef PART_SCAV
1483 & - scav_part*freefelocal +
1484 #else
1485 & - scav*freefelocal +
1486 #endif
1487 & alpfe*inputFelocal/dzlocal
1488 dSidt = - consumpSi + preminSi
1489
1490 c tendency of dissolved organic pool
1491 dDOPdt = totphy_dop + totzoo_dop - DOPremin
1492 dDONdt = totphy_don + totzoo_don - DONremin
1493 dDOFedt = totphy_dofe + totzoo_dofe - DOFeremin
1494 c tendency of particulate detritus pools
1495 dpopdt = totphy_pop + totzoo_pop - preminP + psinkP
1496 dpondt = totphy_pon + totzoo_pon - preminN + psinkN
1497 dpofedt = totphy_pofe + totzoo_pofe - preminFe + psinkFe
1498 dpSidt = totphy_posi + totzoo_posi - preminSi + psinkSi
1499 #ifdef ALLOW_CARBON
1500 dDICdt = - consumpDIC - consumpDIC_PIC
1501 & + preminC + DOCremin
1502 & + disscPIC
1503 dDOCdt = totphy_doc + totzoo_doc - DOCremin
1504 dPOCdt = totphy_poc + totzoo_poc - preminC + psinkC
1505 dPICdt = totphy_pic + totzoo_pic - disscPIC + psinkPIC
1506 dALKdt = - dNO3dt - 2.d0 * (consumpDIC_PIC - disscPIC)
1507 c should be = O2prod - preminP - DOPremin?
1508 c OLD WAY
1509 c dO2dt = - R_OP*dPO4dt
1510 c production of O2 by photosynthesis
1511 dO2dt = R_OP*consumpPO4
1512 c loss of O2 by remineralization
1513 if (O2local.gt.O2crit) then
1514 dO2dt = dO2dt - R_OP*(preminP + DOPremin)
1515 endif
1516 #ifdef OLD_GRAZE
1517 do nz=1,nzmax
1518 dzooCdt(nz) = grazingC(nz)*zooClocal(nz)
1519 & - mortzoo(nz)*zooClocal(nz)
1520 & - mortzoo2(nz)*zooClocal(nz)**2
1521 enddo
1522 #else
1523 do nz=1,nzmax
1524 dzooCdt(nz) = sumgrazzooc(nz)
1525 & - mortzoo(nz)*zooClocal(nz)
1526 & - mortzoo2(nz)*zooClocal(nz)**2
1527 enddo
1528 #endif
1529
1530 #endif
1531
1532 if (debug.eq.10) print*,'dDOPdt', dDOPdt
1533 if (debug.eq.10) print*,'dpopdt',dpopdt
1534 if (debug.eq.10) print*,'dDONdt',dDONdt
1535 if (debug.eq.10) print*,'dpondt',dpondt
1536 c
1537 c -------------------------------------------------------------------
1538 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1539 c --------------------------------------------------------------------------
1540
1541 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-
1542 c Mutation - apply mutation to tendencies [jbmodif]
1543
1544 #ifdef ALLOW_MUTANTS
1545 c apply to all sisters when first sister is encountered
1546 if(runtim .gt. threeyr) then
1547 mutfor=1 _d -8
1548 mutback=1 _d -12
1549 if(numtax .gt. 1)then
1550 do np=1,npro
1551 if(mod(np,numtax).eq. 1. _d 0)then
1552 nsisone = np
1553 nsistwo = np+1
1554 nsisthree = np+2
1555 nsisfour = np+3
1556
1557 grow1 = PspecificPO4(nsisone)
1558 grow2 = PspecificPO4(nsistwo)
1559
1560 if(numtax.eq.2)grow3 = 0.0 _d 0
1561 if(numtax.eq.2)grow4 = 0.0 _d 0
1562
1563 if(numtax.eq.3)grow4 = 0.0 _d 0
1564 if(numtax.ge.3)grow3 = PspecificPO4(nsisthree)
1565
1566 if(numtax.eq.4)grow4 = PspecificPO4(nsisfour)
1567
1568
1569
1570 dphytodt(nsisone) = dphytodt(nsisone)
1571 & - grow1 *1.4427 _d 0*mutfor
1572 & - grow1 *1.4427 _d 0*mutfor
1573 & - grow1 *1.4427 _d 0*mutfor
1574 & + grow2 *1.4427 _d 0*mutback
1575 & + grow3 *1.4427 _d 0*mutback
1576 & + grow4 *1.4427 _d 0*mutback
1577
1578 dphytodt(nsistwo) = dphytodt(nsistwo)
1579 & - grow2 *1.4427 _d 0*mutback
1580 & + grow1 *1.4427 _d 0*mutfor
1581
1582 if(numtax .ge. 3)then
1583 dphytodt(nsisthree) = dphytodt(nsisthree)
1584 & - grow3 *1.4427 _d 0*mutback
1585 & + grow1 *1.4427 _d 0*mutfor
1586 endif
1587
1588 if(numtax .eq. 4)then
1589 dphytodt(nsisfour) = dphytodt(nsisfour)
1590 & - grow4 *1.4427 _d 0*mutback
1591 & + grow1 *1.4427 _d 0*mutfor
1592 c QQQQQQQQQQ FIX FOR NIT RUNS ONLY!!!
1593 if (phyto(nsisfour).eq.0. _d 0) then
1594 if (phyto(nsistwo).eq.0. _d 0) then
1595 if (dphytodt(nsistwo).gt.dphytodt(nsisfour)) then
1596 dphytodt(nsisfour)=dphytodt(nsistwo)
1597 endif
1598 endif
1599 if (phyto(nsisthree).eq.0. _d 0) then
1600 if (dphytodt(nsisthree).gt.dphytodt(nsisfour)) then
1601 dphytodt(nsisfour)=dphytodt(nsisthree)
1602 endif
1603 endif
1604 endif
1605 c QQQQQQQQQQQQQ
1606 endif
1607
1608 c QQQQQQQQQQQQTEST
1609 if (debug.eq.11) then
1610 if (PARlocal.gt.1. _d 0) then
1611 if (dphytodt(nsistwo).gt.dphytodt(nsisfour).and.
1612 & dphytodt(nsisfour).gt.0. _d 0) then
1613 print*,'QQQQ nsistwo>nsisfour',nsistwo,nsisfour,
1614 & dphytodt(nsistwo), dphytodt(nsisfour),
1615 & phyto(nsistwo), phyto(nsisfour),
1616 & phyto(nsisone)
1617 endif
1618 if (dphytodt(nsisthree).gt.dphytodt(nsisfour).and.
1619 & dphytodt(nsisfour).gt.0. _d 0) then
1620 print*,'QQQQ nsisthree>nsisfour',nsisthree,nsisfour,
1621 & dphytodt(nsisthree), dphytodt(nsisfour),
1622 & phyto(nsisthree), phyto(nsisfour),
1623 & phyto(nsisone)
1624 endif
1625 if (dphytodt(nsisfour).gt.dphytodt(nsisone).and.
1626 & dphytodt(nsisone).gt.0. _d 0) then
1627 print*,' BIG QQQQ nsisfour>nsisone',nsisone,nsisfour,
1628 & dphytodt(nsisfour), dphytodt(nsisone),
1629 & phyto(nsisfour), phyto(nsisone)
1630 endif
1631 endif
1632 endif
1633 c QQQQQQQQQTEST
1634 endif
1635 enddo
1636 endif
1637 endif
1638
1639 c mutation is finished
1640 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-
1641 #endif
1642
1643
1644
1645 RETURN
1646 END
1647 #endif /*MONOD*/
1648 #endif /*ALLOW_PTRACERS*/
1649 c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22