/[MITgcm]/MITgcm_contrib/ecco_darwin/v4_3deg/code/darwin_plankton.F
ViewVC logotype

Annotation of /MITgcm_contrib/ecco_darwin/v4_3deg/code/darwin_plankton.F

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


Revision 1.1 - (hide annotations) (download)
Mon Jan 27 02:52:02 2020 UTC (6 years, 7 months ago) by dimitri
Branch: MAIN
First complete check in with v4_llc270 ecco-darwin set-up

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

  ViewVC Help
Powered by ViewVC 1.1.22