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

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

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


Revision 1.1 - (hide annotations) (download)
Tue Jan 14 18:23:29 2020 UTC (5 years, 6 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
- regressing v4_3deg to same vintage code as v4_llc270
- also adding some laggard files in v5_llc270/code_darwin

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

  ViewVC Help
Powered by ViewVC 1.1.22