/[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.6 - (hide annotations) (download)
Thu May 31 21:08:25 2012 UTC (13 years, 1 month ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63n_20120604
Changes since 1.5: +152 -7 lines
o add CDOM-like tracer (#define ALLOW_CDOM)

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

  ViewVC Help
Powered by ViewVC 1.1.22