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

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

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


Revision 1.5 - (hide annotations) (download)
Tue Aug 6 18:28:13 2019 UTC (7 years ago) by dcarroll
Branch: MAIN
Changes since 1.4: +5 -3 lines
Updated V4 code for non-linear dissolution and added additional CO2 flux budget terms

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

  ViewVC Help
Powered by ViewVC 1.1.22