/[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.4 - (hide annotations) (download)
Fri Oct 21 20:56:53 2011 UTC (13 years, 9 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63d_20111107
Changes since 1.3: +30 -12 lines
o add quadratic mortality for zooplankton

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

  ViewVC Help
Powered by ViewVC 1.1.22