/[MITgcm]/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.10 - (hide annotations) (download)
Fri Nov 16 19:50:22 2012 UTC (12 years, 8 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64a_20121116, ctrb_darwin2_ckpt64e_20130305, ctrb_darwin2_ckpt64c_20130120, ctrb_darwin2_ckpt64b_20121224, ctrb_darwin2_ckpt64d_20130219
Changes since 1.9: +43 -21 lines
o Ek, EkoverE diagnostics

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

  ViewVC Help
Powered by ViewVC 1.1.22