/[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.9 - (hide annotations) (download)
Thu Jul 26 18:01:22 2012 UTC (12 years, 11 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt63r_20120817, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt64_20121012
Changes since 1.8: +40 -2 lines
o bug fix, so denitrification (and remineralization) will not happen when
  NO3 below a critical value (NO3crit)

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

  ViewVC Help
Powered by ViewVC 1.1.22