/[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.5 - (hide annotations) (download)
Tue May 19 15:23:46 2015 UTC (10 years, 2 months ago) by benw
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt66m_20171213, HEAD
Changes since 1.4: +19 -35 lines
Ben Ward - modifications to parameters and some structural changes
         - performs much better at OWS Mike

1 benw 1.5 c $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/quota/quota_generate_phyto.F,v 1.4 2015/05/19 14:32:43 benw 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 benw 1.5 ESD1 = 0.5 _d 0 ! minimum plankton ESD
83     ni = 3 ! ni size classes within diameter gaps of * 10
84     factor = 1000. _d 0 ** (1. _d 0 / float(ni))
85 benw 1.4 tau1=1.0 _d 0
86     tau2=1.0 _d 0
87 jahn 1.1 c Allocate Phytoplankton Taxa
88     c Prochloro
89     do jp=1,2
90 benw 1.5 biovol(jp) = pi*ESD1**3/6. _d 0 *factor**(jp-1)
91 jahn 1.1 autotrophy(jp)= 1.00 _d 0
92 benw 1.2 use_NO3(jp) = 1
93 jahn 1.1 use_Si(jp) = 0
94     taxon_mu(jp) = 1.00 _d 0
95     pft(jp) = 1
96     enddo
97     c Synnecho
98 benw 1.5 do jp=3,4
99     biovol(jp) = pi*ESD1**3/6. _d 0 *factor**(jp-2)
100 jahn 1.1 autotrophy(jp)= 1.00 _d 0
101     use_NO3(jp) = 1
102     use_Si(jp) = 0
103 benw 1.2 taxon_mu(jp) = 1.40 _d 0
104 jahn 1.1 pft(jp) = 2
105     enddo
106     c Small Euk
107 benw 1.5 do jp=5,9
108     biovol(jp) = pi*ESD1**3/6. _d 0 *factor**(jp-3)
109 jahn 1.1 autotrophy(jp)= 1.0 _d 0
110     use_NO3(jp) = 1
111     use_Si(jp) = 0
112     taxon_mu(jp) = 2.10 _d 0
113     pft(jp) = 3
114     enddo
115     c Diatoms
116     do jp=10,15
117 benw 1.5 biovol(jp) = pi*ESD1**3/6. _d 0 *factor**(jp-7)
118 jahn 1.1 autotrophy(jp)= 1.0 _d 0
119     use_NO3(jp) = 1
120     use_Si(jp) = 0
121     taxon_mu(jp) = 3.8 _d 0
122     pft(jp) = 4
123     enddo
124     c Specialist grazers
125     do jp=16,16
126 benw 1.5 biovol(jp) = pi*ESD1**3/6. _d 0 *factor**(jp-10)
127 jahn 1.1 autotrophy(jp)= 0.00 _d 0
128     use_NO3(jp) = 0
129     use_Si(jp) = 0
130     taxon_mu(jp) = 0.00 _d 0
131     pft(jp) = 6
132     enddo
133     c
134     do jp=1,16
135     heterotrophy(jp)=1.0 _d 0 - autotrophy(jp)
136     enddo
137 benw 1.4 #else
138     ESD1 = 0.5 _d 0 ! minimum plankton ESD
139     ni = 2 ! ni size classes within diameter gaps of * 10
140     factor = 1000. _d 0 ** (1. _d 0 / float(ni))
141     tau1 = 1.25 _d 0
142     tau2 = 1.0 _d 0 / tau1
143     ntroph = 2
144     if (ntroph.eq.1) then
145     tau1 = 0.0 _d 0
146     tau2 = 0.0 _d 0
147     endif
148     c Allocate plankton traits
149     icount=0
150     do jp2=1,ntroph
151     do jp=1,npmax/ntroph
152     icount = icount + 1
153     biovol(icount) = pi*ESD1**3/6. _d 0 *factor**(jp-1)
154     logvol(icount) = log10(biovol(icount))
155     use_NO3(icount) = 1
156     use_Si(icount) = 0
157     if (ntroph.gt.1) then
158     autotrophy(icount) = 1.0 _d 0 - float(jp2-1)
159     & / float(ntroph - 1)
160     heterotrophy(icount)= 1.0 _d 0 - autotrophy(icount)
161     else
162     autotrophy(icount) = 0.5 _d 0
163     heterotrophy(icount)= 0.5 _d 0
164     endif
165     enddo
166     enddo
167     #endif
168 jahn 1.1 c ----------------------------------------------------------------------
169     c Allometry
170     #ifdef UNCERTAINTY
171     error = 1.0 _d 0
172     #else
173     error = 0.0 _d 0
174     ! set stdev of allometric parameters to zero
175     #endif
176     c ----------------------------------------------------------------------
177     do jp=1,npmax
178     ! parameters independent of nutrient element
179     c CARBON CONTENT
180     p = darwin_random(myThid)
181     call invnormal(a,p,
182     & log10(a_qcarbon),log10(ae_qcarbon)*error)
183     call invnormal(b,p,b_qcarbon,be_qcarbon*error)
184     qcarbon(jp) = 10. _d 0**a * biovol(jp) ** b
185     c INITIAL SLOPE P-I
186     p = darwin_random(myThid)
187     call invnormal(a,p,
188     & log10(a_alphachl),log10(ae_alphachl)*error)
189     call invnormal(b,p,b_alphachl,be_alphachl*error)
190     alphachl(jp) = 10. _d 0**a * biovol(jp) ** b
191     c RESPIRATION RATE
192     p = darwin_random(myThid)
193     IF (a_respir.NE.0. _d 0) THEN
194     call invnormal(a,p,
195     & log10(a_respir),log10(ae_respir)*error)
196     call invnormal(b,p,b_respir,be_respir*error)
197 benw 1.4 respiration(jp) = 10. _d 0**a * biovol(jp) ** b
198 jahn 1.1 ELSE
199     respiration(jp) = 0.0 _d 0
200     ENDIF
201 benw 1.4 c GRAZING SIZE PREFERENCE RATIO
202     p = darwin_random(myThid)
203     call invnormal(a,p,
204     & log10(a_prdpry),log10(ae_prdpry)*error)
205     call invnormal(b,p,b_prdpry,be_prdpry*error)
206     pp_opt(jp) = 10. _d 0**a * biovol(jp) ** b
207     c MAXIMUM GRAZING RATE + WIDTH OF GRAZING KERNEL
208 jahn 1.1 p = darwin_random(myThid)
209     call invnormal(a,p,
210     & log10(a_graz),log10(ae_graz)*error)
211     call invnormal(b,p,b_graz,be_graz*error)
212 benw 1.5 #ifdef ONEGRAZER
213     ! if only one grazer, set max grazing by pp_opt(jp) * prey size
214     graz(jp) = 10. _d 0**a * (biovol(jp)*pp_opt(jp)) ** b
215     #else
216 benw 1.4 ! set grazing rate by grazer size, non-grazers to zero
217     graz(jp) = 10. _d 0**a * biovol(jp) ** b
218     & * heterotrophy(jp) ** tau2
219 benw 1.5 #endif
220     pp_sig(jp) = 2. _d 0
221     c FRACTION GRAZED AND MORTALITY TO DOM
222 jahn 1.1 do io=1,iomax-iChl
223 benw 1.4 #ifdef ALLOWPFT
224     if (pft(jp).lt.3) beta_graz(io,jp)=0.8
225     if (pft(jp).gt.2) beta_graz(io,jp)=0.5
226     #else
227     beta_graz(io,jp) = 0.9 _d 0 - 0.7 _d 0
228     & / (1.0 _d 0 + exp(-logvol(jp)+2.0 _d 0))
229     #endif
230     beta_mort(io,jp) = beta_graz(io,jp)
231 jahn 1.1 enddo
232     c GRAZING HALF-SATURATION
233     p = darwin_random(myThid)
234     call invnormal(a,p,
235     & log10(a_kg),log10(ae_kg)*error)
236     call invnormal(b,p,b_kg,be_kg*error)
237     kg(jp) = 10. _d 0**a * biovol(jp) ** b
238 benw 1.4 #ifdef DIFFLIMIT
239     & * heterotrophy(jp) ** tau1
240     #endif
241 jahn 1.1 c PHYTOPLANKTON SINKING
242     p = darwin_random(myThid)
243     call invnormal(a,p,
244     & log10(a_biosink),log10(ae_biosink)*error)
245     call invnormal(b,p,b_biosink,be_biosink*error)
246 benw 1.4 biosink(jp) = (10.0 _d 0**a) * biovol(jp) ** b
247     #ifdef ALLOWPFT
248     if (pft(jp).eq.6) biosink(jp) = 0. _d 0
249     #endif
250     c MORTALITY
251     ! constant background mortality
252 jahn 1.1 p = darwin_random(myThid)
253     call invnormal(a,p,
254 benw 1.4 & log10(a_mort),log10(ae_mort)*error)
255     call invnormal(b,p,b_mort,be_mort*error)
256     kmort(jp) = (10.0 _d 0**a) * biovol(jp) ** b
257 jahn 1.1 ! parameters relating to inorganic nutrients
258     do ii=1,iimax
259     c MAXIMUM NUTRIENT UPTAKE RATE
260     p = darwin_random(myThid)
261     call invnormal(a,p,
262     & log10(a_vmaxi(ii)),log10(ae_vmaxi(ii))*error)
263     call invnormal(b,p,b_vmaxi(ii),be_vmaxi(ii)*error)
264     if (ii.eq.iDIC) then
265 benw 1.4 #ifdef ALLOWPFT
266 jahn 1.1 vmaxi(ii,jp)= 10. _d 0**a * biovol(jp) ** b
267     & * taxon_mu(jp)
268 benw 1.4 #else
269     vmaxi(ii,jp)= (3.1 _d 0+logvol(jp))
270     & / (5.0 _d 0 - 3.8 _d 0*logvol(jp) + logvol(jp)**2)
271     & / 86400. _d 0
272     & * autotrophy(jp) ** tau1
273     #endif
274 jahn 1.1 else
275     vmaxi(ii,jp)= 10. _d 0**a * biovol(jp) ** b
276 benw 1.4 & * autotrophy(jp) ** tau1
277     c NUTRIENT HALF-SATURATION CONSTANT
278 jahn 1.1 p = darwin_random(myThid)
279     call invnormal(a,p,
280     & log10(a_kn(ii)),log10(ae_kn(ii))*error)
281     call invnormal(b,p,b_kn(ii),be_kn(ii)*error)
282     kn(ii,jp) = 10. _d 0**a * biovol(jp) ** b
283 benw 1.4 #ifdef DIFFLIMIT
284     & * autotrophy(jp) ** tau1
285     #endif
286 jahn 1.1 endif
287     enddo
288     #ifdef SQUOTA
289     ! Silicate parameters to zero for non-diatoms
290     vmaxi(iSi,jp) = vmaxi(iSi,jp) * float(use_Si(jp))
291     #endif
292     c
293     if (use_NO3(jp).eq.0) then
294     ! prochlorocococcus can't use NO3
295     vmaxi(iNO3,jp) = 0.0 _d 0
296     ! but have higher NH4 affinity
297     vmaxi(iNH4,jp) = vmaxi(iNH4,jp) * 2.0 _d 0
298     endif
299     ! parameters relating to quota nutrients
300     do io=1,iomax-iChl
301     c EXCRETION
302     if ((io.eq.iCarb.or.io.eq.iNitr.or.io.eq.iPhos)
303     & .and.a_kexc(io).NE.0. _d 0
304     & .and.ae_kexc(io).NE.0. _d 0) then
305     p = darwin_random(myThid)
306     call invnormal(a,p,
307     & log10(a_kexc(io)),log10(ae_kexc(io))*error)
308     call invnormal(b,p,b_kexc(io),be_kexc(io)*error)
309     kexc(io,jp) = 10. _d 0**a * biovol(jp) ** b
310     else
311     kexc(io,jp) = 0. _d 0
312     endif
313     if (io.ne.iCarb) then
314     c MINIMUM QUOTA
315     p = darwin_random(myThid)
316     call invnormal(a,p,
317     & log10(a_qmin(io)),log10(ae_qmin(io))*error)
318     call invnormal(b,p,b_qmin(io),be_qmin(io)*error)
319     qmin(io,jp) = 10. _d 0**a * biovol(jp) ** b
320 benw 1.4 ! & * (autotrophy(jp) ** tau2
321     ! & + heterotrophy(jp) ** tau2)
322 jahn 1.1 c MAXIMUM QUOTA
323     p = darwin_random(myThid)
324     call invnormal(a,p,
325     & log10(a_qmax(io)),log10(ae_qmax(io))*error)
326     call invnormal(b,p,b_qmax(io),be_qmax(io)*error)
327     qmax(io,jp) = 10. _d 0**a * biovol(jp) ** b
328     endif
329     enddo
330     #ifdef SQUOTA
331     ! Silicate parameters to zero for non-diatoms
332     qmin(iSili,jp) = qmin(iSili,jp) * float(use_Si(jp))
333     qmax(iSili,jp) = qmax(iSili,jp) * float(use_Si(jp))
334     #endif
335     c PREFERENCE FUNCTION
336     ! assign grazing preference according to predator/prey radius ratio
337     do jp2=1,npmax ! jp2 denotes prey
338 benw 1.4 if (heterotrophy(jp).gt.0. _d 0) then
339 jahn 1.1 prd_pry = biovol(jp) / biovol(jp2)
340     graz_pref(jp,jp2) =
341 benw 1.5 #ifdef ONEGRAZER
342 benw 1.2 & 1.0 _d 0
343 jahn 1.1 #else
344 benw 1.4 & exp(-log(prd_pry/pp_opt(jp))**2 / (2*pp_sig(jp)**2))
345 jahn 1.1 #endif
346     if (graz_pref(jp,jp2).lt.1. _d -4) then
347     graz_pref(jp,jp2)=0. _d 0
348     endif
349     assim_graz(jp,jp2) = ass_eff
350     else
351     graz_pref(jp,jp2) = 0. _d 0
352     endif
353     enddo
354     c
355     c..........................................................
356     c generate phyto Temperature Function parameters
357     c.......................................................
358     phytoTempCoeff(jp) = tempcoeff1
359     phytoTempExp1(jp) = tempcoeff3
360     phytoTempExp2(jp) = tempcoeff2_small
361     & + (tempcoeff2_big-tempcoeff2_small)
362     & * float(jp-1)/npmaxm1
363     phytoTempOptimum(jp) = 2. _d 0
364     phytoDecayPower(jp) = tempdecay
365    
366     c..........................................................
367     enddo
368    
369    
370     RETURN
371     END
372     #endif /*ALLOW_QUOTA*/
373     #endif /*ALLOW_DARWIN*/
374     #endif /*ALLOW_PTRACERS*/
375    
376     c ===========================================================

  ViewVC Help
Powered by ViewVC 1.1.22