/[MITgcm]/MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/monod/monod_plankton.F

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


Revision 1.5 - (hide annotations) (download)
Tue Nov 15 20:16:34 2011 UTC (13 years, 8 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63k_20120317
Changes since 1.4: +1 -2 lines
remove unnused header

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

  ViewVC Help
Powered by ViewVC 1.1.22