/[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.15 - (hide annotations) (download)
Thu May 1 16:19:32 2014 UTC (11 years, 2 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt65j_20150225, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt65_20140718, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt65e_20140929, ctrb_darwin2_ckpt64z_20140711, ctrb_darwin2_ckpt65l_20150504, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt64y_20140622, ctrb_darwin2_ckpt65d_20140915, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt64x_20140524, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt65g_20141120, ctrb_darwin2_ckpt65k_20150402, ctrb_darwin2_ckpt64w_20140502, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt65f_20141014, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt65i_20150123, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt65a_20140728, ctrb_darwin2_ckpt65b_20140812, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt65c_20140830, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt65h_20141217, ctrb_darwin2_ckpt66m_20171213, HEAD
Changes since 1.14: +35 -1 lines
add option FIX_ZOO_QUOTAS

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

  ViewVC Help
Powered by ViewVC 1.1.22