/[MITgcm]/MITgcm_contrib/darwin2/pkg/quota/quota_generate_phyto.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/quota/quota_generate_phyto.F

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


Revision 1.4 - (hide annotations) (download)
Tue May 19 14:32:43 2015 UTC (10 years, 2 months ago) by benw
Branch: MAIN
Changes since 1.3: +111 -86 lines
Ben Ward - some superficial structural changes allowing runs with no pfts
         - more significant structural and parameter changes to follow later

1 benw 1.4 c $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/quota/quota_generate_phyto.F,v 1.3 2013/06/12 17:53:27 jahn Exp $
2 jahn 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5     #include "PTRACERS_OPTIONS.h"
6     #include "DARWIN_OPTIONS.h"
7    
8     #ifdef ALLOW_PTRACERS
9     #ifdef ALLOW_DARWIN
10     #ifdef ALLOW_QUOTA
11    
12     c ==========================================================
13     c SUBROUTINE QUOTA_GENERATE_PHYTO
14     c generate parameters for "Operational Taxonomic Units" of plankton (index jp)
15     c using an allometric approach
16     c
17     c Ben Ward 2009/10
18     c ==========================================================
19     SUBROUTINE QUOTA_GENERATE_PHYTO(myThid)
20    
21     implicit none
22     #include "EEPARAMS.h"
23     #include "DARWIN_PARAMS.h"
24     #include "QUOTA_SIZE.h"
25     #include "QUOTA.h"
26    
27     C !INPUT PARAMETERS: ===================================================
28     C myThid :: thread number
29     INTEGER myThid
30    
31     C === Functions ===
32     _RL DARWIN_RANDOM
33     EXTERNAL DARWIN_RANDOM
34     _RL DARWIN_RANDOM_NORMAL
35     EXTERNAL DARWIN_RANDOM_NORMAL
36    
37    
38     C !LOCAL VARIABLES:
39     C === Local variables ===
40     C msgBuf - Informational/error meesage buffer
41     CHARACTER*(MAX_LEN_MBUF) msgBuf
42    
43     _RL RandNo
44     _RL mortdays
45     _RL year
46     _RL rtime
47     _RL standin
48     _RL tmpsrt
49     _RL tmpend
50     _RL tmprng
51     _RL iimaxm1
52     _RL npmaxm1
53     _RL komaxm1
54     _RL prd_pry
55     _RL factor
56 benw 1.4 #ifdef ALLOWPFT
57 jahn 1.1 _RL taxon_mu(npmax)
58 benw 1.4 #endif
59 jahn 1.1 _RL a,b,p,error
60     _RL heterotrophy(npmax)
61 benw 1.4 _RL tau1,tau2
62     _RL ESD1,pi
63     _RL logvol(npmax)
64     INTEGER ii,io,jp,ko,ni
65     INTEGER jp2,icount,ntroph
66 jahn 1.1 INTEGER signvar
67     CEOP
68     c
69 benw 1.4 standin= 0. _d 0
70     pi = 4. _d 0 * datan(1. _d 0)
71 jahn 1.1 c each time generate another functional group add one to ngroups
72     ngroups = ngroups + 1
73    
74     iimaxm1 = float(iimax-1)
75     npmaxm1 = float(npmax-1)
76     komaxm1 = float(komax-1)
77     c
78     c..........................................................
79     c Generate plankton volumes and stochastic parameters
80     c..........................................................
81 benw 1.4 #ifdef ALLOWPFT
82 jahn 1.1 factor = 2. _d 0
83 benw 1.4 tau1=1.0 _d 0
84     tau2=1.0 _d 0
85 jahn 1.1 c Allocate Phytoplankton Taxa
86     c Prochloro
87     do jp=1,2
88 benw 1.2 biovol(jp) = 0.125 _d 0 * factor**(jp-1)
89 jahn 1.1 autotrophy(jp)= 1.00 _d 0
90 benw 1.2 use_NO3(jp) = 1
91 jahn 1.1 use_Si(jp) = 0
92     taxon_mu(jp) = 1.00 _d 0
93     pft(jp) = 1
94     enddo
95     c Synnecho
96     do jp=3,5
97     biovol(jp) = 0.50 _d 0 * factor**(jp-3)
98     autotrophy(jp)= 1.00 _d 0
99     use_NO3(jp) = 1
100     use_Si(jp) = 0
101 benw 1.2 taxon_mu(jp) = 1.40 _d 0
102 jahn 1.1 pft(jp) = 2
103     enddo
104     c Small Euk
105     do jp=6,9
106     biovol(jp) = 4.00 _d 0 * factor**(jp-6)
107     autotrophy(jp)= 1.0 _d 0
108     use_NO3(jp) = 1
109     use_Si(jp) = 0
110     taxon_mu(jp) = 2.10 _d 0
111     pft(jp) = 3
112     enddo
113     c Diatoms
114     do jp=10,15
115 jahn 1.3 biovol(jp) = 128.0 _d 0 * factor**(jp-10)
116 jahn 1.1 autotrophy(jp)= 1.0 _d 0
117     use_NO3(jp) = 1
118     use_Si(jp) = 0
119     taxon_mu(jp) = 3.8 _d 0
120     pft(jp) = 4
121     enddo
122     c Specialist grazers
123     do jp=16,16
124 jahn 1.3 biovol(jp) = 8.0 _d 0 * factor**(jp-16)
125 jahn 1.1 autotrophy(jp)= 0.00 _d 0
126     use_NO3(jp) = 0
127     use_Si(jp) = 0
128     taxon_mu(jp) = 0.00 _d 0
129     pft(jp) = 6
130     enddo
131     c
132     do jp=1,16
133     heterotrophy(jp)=1.0 _d 0 - autotrophy(jp)
134     enddo
135 benw 1.4 #else
136     ESD1 = 0.5 _d 0 ! minimum plankton ESD
137     ni = 2 ! ni size classes within diameter gaps of * 10
138     factor = 1000. _d 0 ** (1. _d 0 / float(ni))
139     tau1 = 1.25 _d 0
140     tau2 = 1.0 _d 0 / tau1
141     ntroph = 2
142     if (ntroph.eq.1) then
143     tau1 = 0.0 _d 0
144     tau2 = 0.0 _d 0
145     endif
146     c Allocate plankton traits
147     icount=0
148     do jp2=1,ntroph
149     do jp=1,npmax/ntroph
150     icount = icount + 1
151     biovol(icount) = pi*ESD1**3/6. _d 0 *factor**(jp-1)
152     logvol(icount) = log10(biovol(icount))
153     use_NO3(icount) = 1
154     use_Si(icount) = 0
155     if (ntroph.gt.1) then
156     autotrophy(icount) = 1.0 _d 0 - float(jp2-1)
157     & / float(ntroph - 1)
158     heterotrophy(icount)= 1.0 _d 0 - autotrophy(icount)
159     else
160     autotrophy(icount) = 0.5 _d 0
161     heterotrophy(icount)= 0.5 _d 0
162     endif
163     enddo
164     enddo
165     #endif
166 jahn 1.1 c ----------------------------------------------------------------------
167     c Allometry
168     #ifdef UNCERTAINTY
169     error = 1.0 _d 0
170     #else
171     error = 0.0 _d 0
172     ! set stdev of allometric parameters to zero
173     #endif
174     c ----------------------------------------------------------------------
175     do jp=1,npmax
176     ! parameters independent of nutrient element
177     c CARBON CONTENT
178     p = darwin_random(myThid)
179     call invnormal(a,p,
180     & log10(a_qcarbon),log10(ae_qcarbon)*error)
181     call invnormal(b,p,b_qcarbon,be_qcarbon*error)
182     qcarbon(jp) = 10. _d 0**a * biovol(jp) ** b
183     c INITIAL SLOPE P-I
184     p = darwin_random(myThid)
185     call invnormal(a,p,
186     & log10(a_alphachl),log10(ae_alphachl)*error)
187     call invnormal(b,p,b_alphachl,be_alphachl*error)
188     alphachl(jp) = 10. _d 0**a * biovol(jp) ** b
189     c RESPIRATION RATE
190     p = darwin_random(myThid)
191     IF (a_respir.NE.0. _d 0) THEN
192     call invnormal(a,p,
193     & log10(a_respir),log10(ae_respir)*error)
194     call invnormal(b,p,b_respir,be_respir*error)
195 benw 1.4 respiration(jp) = 10. _d 0**a * biovol(jp) ** b
196 jahn 1.1 ELSE
197     respiration(jp) = 0.0 _d 0
198     ENDIF
199 benw 1.4 c GRAZING SIZE PREFERENCE RATIO
200     p = darwin_random(myThid)
201     call invnormal(a,p,
202     & log10(a_prdpry),log10(ae_prdpry)*error)
203     call invnormal(b,p,b_prdpry,be_prdpry*error)
204     pp_opt(jp) = 10. _d 0**a * biovol(jp) ** b
205     c MAXIMUM GRAZING RATE + WIDTH OF GRAZING KERNEL
206 jahn 1.1 p = darwin_random(myThid)
207     call invnormal(a,p,
208     & log10(a_graz),log10(ae_graz)*error)
209     call invnormal(b,p,b_graz,be_graz*error)
210 benw 1.4 !#ifdef ONEGRAZER
211     !! if only one grazer, set max grazing by 1 * prey size
212     !! only P are grazed
213     ! graz(jp) = 10. _d 0**a * (biovol(jp)*pp_opt(jp)) ** b
214     ! & * autotrophy(jp) ** tau2
215     !! temp fix to remove size dep.
216     ! graz(jp) = 10. _d 0**a * (biovol(npmax)) ** b
217     ! & * autotrophy(jp) ** tau2
218     !! temp fix to remove size dep.
219     ! pp_sig(jp) = 1. _d 10
220     !#else
221     ! set grazing rate by grazer size, non-grazers to zero
222     graz(jp) = 10. _d 0**a * biovol(jp) ** b
223     & * heterotrophy(jp) ** tau2
224     pp_sig(jp) = 2. _d 0
225     !#endif
226     c FRACTION GRAZED TAND MORTALITY TO DOM
227 jahn 1.1 do io=1,iomax-iChl
228 benw 1.4 #ifdef ALLOWPFT
229     if (pft(jp).lt.3) beta_graz(io,jp)=0.8
230     if (pft(jp).gt.2) beta_graz(io,jp)=0.5
231     #else
232     beta_graz(io,jp) = 0.9 _d 0 - 0.7 _d 0
233     & / (1.0 _d 0 + exp(-logvol(jp)+2.0 _d 0))
234     #endif
235     beta_mort(io,jp) = beta_graz(io,jp)
236 jahn 1.1 enddo
237     c GRAZING HALF-SATURATION
238     p = darwin_random(myThid)
239     call invnormal(a,p,
240     & log10(a_kg),log10(ae_kg)*error)
241     call invnormal(b,p,b_kg,be_kg*error)
242     kg(jp) = 10. _d 0**a * biovol(jp) ** b
243 benw 1.4 #ifdef DIFFLIMIT
244     & * heterotrophy(jp) ** tau1
245     #endif
246 jahn 1.1 c PHYTOPLANKTON SINKING
247     p = darwin_random(myThid)
248     call invnormal(a,p,
249     & log10(a_biosink),log10(ae_biosink)*error)
250     call invnormal(b,p,b_biosink,be_biosink*error)
251 benw 1.4 biosink(jp) = (10.0 _d 0**a) * biovol(jp) ** b
252     #ifdef ALLOWPFT
253     if (pft(jp).eq.6) biosink(jp) = 0. _d 0
254     #endif
255     c MORTALITY
256     ! constant background mortality
257 jahn 1.1 p = darwin_random(myThid)
258     call invnormal(a,p,
259 benw 1.4 & log10(a_mort),log10(ae_mort)*error)
260     call invnormal(b,p,b_mort,be_mort*error)
261     kmort(jp) = (10.0 _d 0**a) * biovol(jp) ** b
262 jahn 1.1 ! parameters relating to inorganic nutrients
263     do ii=1,iimax
264     c MAXIMUM NUTRIENT UPTAKE RATE
265     p = darwin_random(myThid)
266     call invnormal(a,p,
267     & log10(a_vmaxi(ii)),log10(ae_vmaxi(ii))*error)
268     call invnormal(b,p,b_vmaxi(ii),be_vmaxi(ii)*error)
269     if (ii.eq.iDIC) then
270 benw 1.4 #ifdef ALLOWPFT
271 jahn 1.1 vmaxi(ii,jp)= 10. _d 0**a * biovol(jp) ** b
272     & * taxon_mu(jp)
273 benw 1.4 #else
274     vmaxi(ii,jp)= (3.1 _d 0+logvol(jp))
275     & / (5.0 _d 0 - 3.8 _d 0*logvol(jp) + logvol(jp)**2)
276     & / 86400. _d 0
277     & * autotrophy(jp) ** tau1
278     #endif
279 jahn 1.1 else
280     vmaxi(ii,jp)= 10. _d 0**a * biovol(jp) ** b
281 benw 1.4 & * autotrophy(jp) ** tau1
282     c NUTRIENT HALF-SATURATION CONSTANT
283 jahn 1.1 p = darwin_random(myThid)
284     call invnormal(a,p,
285     & log10(a_kn(ii)),log10(ae_kn(ii))*error)
286     call invnormal(b,p,b_kn(ii),be_kn(ii)*error)
287     kn(ii,jp) = 10. _d 0**a * biovol(jp) ** b
288 benw 1.4 #ifdef DIFFLIMIT
289     & * autotrophy(jp) ** tau1
290     #endif
291 jahn 1.1 endif
292     enddo
293     #ifdef SQUOTA
294     ! Silicate parameters to zero for non-diatoms
295     vmaxi(iSi,jp) = vmaxi(iSi,jp) * float(use_Si(jp))
296     #endif
297     c
298     if (use_NO3(jp).eq.0) then
299     ! prochlorocococcus can't use NO3
300     vmaxi(iNO3,jp) = 0.0 _d 0
301     ! but have higher NH4 affinity
302     vmaxi(iNH4,jp) = vmaxi(iNH4,jp) * 2.0 _d 0
303     endif
304     ! parameters relating to quota nutrients
305     do io=1,iomax-iChl
306     c EXCRETION
307     if ((io.eq.iCarb.or.io.eq.iNitr.or.io.eq.iPhos)
308     & .and.a_kexc(io).NE.0. _d 0
309     & .and.ae_kexc(io).NE.0. _d 0) then
310     p = darwin_random(myThid)
311     call invnormal(a,p,
312     & log10(a_kexc(io)),log10(ae_kexc(io))*error)
313     call invnormal(b,p,b_kexc(io),be_kexc(io)*error)
314     kexc(io,jp) = 10. _d 0**a * biovol(jp) ** b
315     else
316     kexc(io,jp) = 0. _d 0
317     endif
318     if (io.ne.iCarb) then
319     c MINIMUM QUOTA
320     p = darwin_random(myThid)
321     call invnormal(a,p,
322     & log10(a_qmin(io)),log10(ae_qmin(io))*error)
323     call invnormal(b,p,b_qmin(io),be_qmin(io)*error)
324     qmin(io,jp) = 10. _d 0**a * biovol(jp) ** b
325 benw 1.4 ! & * (autotrophy(jp) ** tau2
326     ! & + heterotrophy(jp) ** tau2)
327 jahn 1.1 c MAXIMUM QUOTA
328     p = darwin_random(myThid)
329     call invnormal(a,p,
330     & log10(a_qmax(io)),log10(ae_qmax(io))*error)
331     call invnormal(b,p,b_qmax(io),be_qmax(io)*error)
332     qmax(io,jp) = 10. _d 0**a * biovol(jp) ** b
333     endif
334     enddo
335     #ifdef SQUOTA
336     ! Silicate parameters to zero for non-diatoms
337     qmin(iSili,jp) = qmin(iSili,jp) * float(use_Si(jp))
338     qmax(iSili,jp) = qmax(iSili,jp) * float(use_Si(jp))
339     #endif
340 benw 1.4 #ifdef ALLOWPFT
341 jahn 1.1 c Zooplankton have approximately Redfieldian N:C ratio
342     if (pft(jp).eq.6) then
343     qmin(iNitr,jp) = 0.0755 _d 0
344     qmax(iNitr,jp) = 0.1510 _d 0
345     endif
346 benw 1.4 #endif
347 jahn 1.1 c PREFERENCE FUNCTION
348     ! assign grazing preference according to predator/prey radius ratio
349     do jp2=1,npmax ! jp2 denotes prey
350 benw 1.4 #ifdef ALLOWPFT
351 jahn 1.1 if (heterotrophy(jp).gt.0. _d 0.and.pft(jp2).ne.6) then
352 benw 1.4 #else
353     if (heterotrophy(jp).gt.0. _d 0) then
354     #endif
355 jahn 1.1 prd_pry = biovol(jp) / biovol(jp2)
356     graz_pref(jp,jp2) =
357 benw 1.4 #ifdef SWITCHING
358 benw 1.2 & 1.0 _d 0
359 jahn 1.1 #else
360 benw 1.4 & exp(-log(prd_pry/pp_opt(jp))**2 / (2*pp_sig(jp)**2))
361 jahn 1.1 #endif
362     if (graz_pref(jp,jp2).lt.1. _d -4) then
363     graz_pref(jp,jp2)=0. _d 0
364     endif
365     assim_graz(jp,jp2) = ass_eff
366     else
367     graz_pref(jp,jp2) = 0. _d 0
368     endif
369     enddo
370     c
371     c..........................................................
372     c generate phyto Temperature Function parameters
373     c.......................................................
374     phytoTempCoeff(jp) = tempcoeff1
375     phytoTempExp1(jp) = tempcoeff3
376     phytoTempExp2(jp) = tempcoeff2_small
377     & + (tempcoeff2_big-tempcoeff2_small)
378     & * float(jp-1)/npmaxm1
379     phytoTempOptimum(jp) = 2. _d 0
380     phytoDecayPower(jp) = tempdecay
381    
382     c..........................................................
383     enddo
384    
385    
386     RETURN
387     END
388     #endif /*ALLOW_QUOTA*/
389     #endif /*ALLOW_DARWIN*/
390     #endif /*ALLOW_PTRACERS*/
391    
392     c ===========================================================

  ViewVC Help
Powered by ViewVC 1.1.22