/[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.7 - (hide annotations) (download)
Thu Jun 28 20:36:11 2012 UTC (13 years ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63o_20120629
Changes since 1.6: +98 -46 lines
o add temperature dependencies for mortality (default to no dependence - need to
  change comments in monod_tempfunc to enable)

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

  ViewVC Help
Powered by ViewVC 1.1.22