/[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.1 - (hide annotations) (download)
Thu Dec 21 22:16:32 2017 UTC (8 years, 8 months ago) by dimitri
Branch: MAIN
Making the bottom boundary porous to particulates.
In other words, particulates are removed from
suspension and assumed to accumulate in sediments.

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

  ViewVC Help
Powered by ViewVC 1.1.22