/[MITgcm]/MITgcm_contrib/ecco_darwin/v4_llc270/code_darwin/darwin_plankton.F
ViewVC logotype

Annotation of /MITgcm_contrib/ecco_darwin/v4_llc270/code_darwin/darwin_plankton.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Wed Jun 5 02:36:33 2019 UTC (7 years, 3 months ago) by dimitri
Branch: MAIN
Changes since 1.1: +19 -4 lines
Adding Naviaux (2019, Marine Chemistry) dissolution rate.

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

  ViewVC Help
Powered by ViewVC 1.1.22