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

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

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


Revision 1.1 - (hide annotations) (download)
Wed Apr 13 18:56:25 2011 UTC (14 years, 3 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_baseline
darwin2 initial checkin

1 jahn 1.1 C $Header$
2     C $Name$
3    
4     #include "CPP_OPTIONS.h"
5     #include "DARWIN_OPTIONS.h"
6    
7     #ifdef ALLOW_PTRACERS
8     #ifdef ALLOW_MONOD
9    
10     c ==========================================================
11     c SUBROUTINE MONOD_INIT_VARI()
12     c initialize stuff for generalized plankton model
13     c adapted from NPZD2Fe - Mick Follows, Fall 2005
14     c modified - Stephanie Dutkiewicz, Spring 2006
15     c ==========================================================
16     c
17     SUBROUTINE MONOD_INIT_VARI(myThid)
18    
19     IMPLICIT NONE
20    
21     #include "SIZE.h"
22     #include "GRID.h"
23     #include "DYNVARS.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "MONOD_SIZE.h"
27     #include "MONOD.h"
28     #include "DARWIN_IO.h"
29    
30     c ANNA define params for WAVEBANDS
31     #ifdef WAVEBANDS
32     #include "SPECTRAL_SIZE.h"
33     #include "SPECTRAL.h"
34     #include "WAVEBANDS_PARAMS.h"
35     #endif
36    
37     #ifdef ALLOW_DIAZ
38     #include "PTRACERS_SIZE.h"
39     #include "PTRACERS_PARAMS.h"
40     #include "PTRACERS_FIELDS.h"
41     #endif
42    
43    
44     C !INPUT PARAMETERS: ===================================================
45     C myThid :: thread number
46     INTEGER myThid
47    
48     C === Functions ===
49     _RL DARWIN_RANDOM
50     EXTERNAL DARWIN_RANDOM
51    
52     C !LOCAL VARIABLES:
53     C === Local variables ===
54     C msgBuf - Informational/error meesage buffer
55     CHARACTER*(MAX_LEN_MBUF) msgBuf
56     INTEGER IniUnit1, IniUnit2
57    
58     INTEGER bi, bj, k, i, j, iPAR
59    
60     INTEGER np
61     INTEGER nz
62     c ANNA need nl for wavebands
63     #ifdef WAVEBANDS
64     integer ilam
65     integer nl
66     _RL cu_area
67     #endif
68    
69     CEOP
70    
71     WRITE(msgBuf,'(A)')
72     & '// ======================================================='
73     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
74     & SQUEEZE_RIGHT, myThid )
75     WRITE(msgBuf,'(A)') '// Darwin init variables >>> START <<<'
76     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
77     & SQUEEZE_RIGHT, myThid )
78     WRITE(msgBuf,'(A)')
79     & '// ======================================================='
80     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
81     & SQUEEZE_RIGHT, myThid )
82    
83     c test....................
84     c write(6,*)'testing in npzd2fe_init_vari '
85     c test....................
86    
87    
88     c set up ecosystem coefficients
89     c
90     c seed randomization
91     CALL DARWIN_RANDOM_INIT(darwin_seed, myThid)
92    
93     c initialize total number of functional groups tried
94     ngroups = 0
95     do np = 1, npmax
96     #ifdef ALLOW_MUTANTS
97     call MONOD_GENERATE_MUTANTS(MyThid, np)
98     #else
99     call MONOD_GENERATE_PHYTO(MyThid, np)
100     #endif
101     end do
102    
103     c reduce amount of diaz
104     #ifdef ALLOW_DIAZ
105     do np = 1, npmax
106     if (diazotroph(np) .eq. 1. _d 0) then
107     DO bj = myByLo(myThid), myByHi(myThid)
108     DO bi = myBxLo(myThid), myBxHi(myThid)
109     DO j=1-Oly,sNy+Oly
110     DO i=1-Olx,sNx+Olx
111     DO k=1,nR
112     Ptracer(i,j,k,bi,bj,iPhy+np-1) =
113     & Ptracer(i,j,k,bi,bj,iPhy+np-1)/10. _d 0
114     ENDDO
115     ENDDO
116     ENDDO
117     ENDDO
118     ENDDO
119     endif
120     enddo
121     #endif
122    
123     c initialize zooplankton
124     call MONOD_GENERATE_ZOO(MyThid)
125    
126     c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127    
128     c ANNA call WAVEBANDS_INIT_VARI to assign variable parameters
129     #ifdef WAVEBANDS
130     call WAVEBANDS_INIT_VARI(MyThid)
131     #endif
132    
133     c ANNA get alphachl from mQyield / aphy_chl
134     c ANNA must do this after params are assigned, but before written out
135     c ANNA use aphy_chl_ps for growth. To turn off, simply set same coefs as aphy_chl in input files.
136     #ifdef GEIDER
137     #ifndef WAVEBANDS
138     do np = 1,npmax
139     alphachl(np) = mQyield(np) * aphy_chl_ave
140     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
141     c no temp/nutrient limitation
142     c chl2cmin(np)=chl2cmax(np)/
143     c & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
144     c & (2*pcmax(np)))
145     chl2cmin(np)=0. _d 0
146     enddo
147     #else
148     do np = 1,npmax
149     do nl = 1,tlam
150     alphachl_nl(np,nl) = mQyield(np) * aphy_chl_ps(np,nl)
151     end do
152     c find mean
153     cu_area = 0.d0
154     do nl = 1,tlam
155     cu_area = cu_area + wb_width(nl)*alphachl_nl(np,nl)
156     end do
157     alpha_mean(np) = cu_area / wb_totalWidth
158    
159     chl2cmin(np)=chl2cmax(np)/
160     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
161     & (2*pcmax(np)))
162     end do
163     #endif
164     #ifdef DYNAMIC_CHL
165     c check Chl fields are reasonable
166     do np=1,npmax
167     chl2cmin(np)=chl2cmax(np)/
168     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
169     & (2*pcmax(np)))
170     enddo
171     call MONOD_CHECK_CHL(myThid)
172     #endif
173     #endif
174     c ANNA endif
175    
176     _BEGIN_MASTER(myThid)
177     c write out initial phyto characteristics
178     #ifndef GEIDER
179     CALL MDSFINDUNIT( IniUnit1, mythid )
180     open(IniUnit1,file='plankton-ini-char.dat',status='unknown')
181     CALL MDSFINDUNIT( IniUnit2, mythid )
182     open(IniUnit2,file='plankton_ini_char_nohead.dat',
183     & status='unknown')
184     #ifdef OLD_GRAZE
185     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip
186     & wsink KsP KsN KsFe KsSi g1 g2 Kpar Kinh
187     & Topt nsrc np'
188     do np = 1, npmax
189     write(IniUnit1,110)diacoc(np),diazotroph(np),physize(np),
190     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
191     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
192     & wsink(np),
193     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
194     & ,KsatSi(np),
195     & graze(np,1),graze(np,2),
196     & KsatPAR(np),Kinhib(np),
197     & phytoTempOptimum(np),nsource(np),np
198     write(IniUnit2,110)diacoc(np),diazotroph(np),physize(np),
199     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
200     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
201     & wsink(np),
202     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
203     & ,KsatSi(np),
204     & graze(np,1),graze(np,2),
205     & KsatPAR(np),Kinhib(np),
206     & phytoTempOptimum(np),nsource(np),np
207     end do
208     #else
209     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip wsink
210     & KsP KsN KsFe KsSi palat1 palat2 Kpar Kinh Topt nsrc
211     & np'
212     do np = 1, npmax
213     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
214     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
215     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
216     & wsink(np),
217     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
218     & ,KsatSi(np),
219     & palat(np,1),palat(np,2),
220     & KsatPAR(np),Kinhib(np),
221     & phytoTempOptimum(np),nsource(np),np
222     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
223     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
224     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
225     & wsink(np),
226     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
227     & ,KsatSi(np),
228     & palat(np,1),palat(np,2),
229     & KsatPAR(np),Kinhib(np),
230     & phytoTempOptimum(np),nsource(np),np
231     end do
232     #endif
233     #endif
234    
235     #ifdef GEIDER
236     c ANNA outputs mQyield as 10^(4) mmol C (uEin)-1
237     CALL MDSFINDUNIT( IniUnit1, mythid )
238     open(IniUnit1,file='gplankton-ini-char.dat',status='unknown')
239     CALL MDSFINDUNIT( IniUnit2, mythid )
240     open(IniUnit2,file='gplankton_ini_char_nohead.dat',
241     & status='unknown')
242     write(IniUnit1,*)'dico diaz size pcmax mort Rnp Rfep Rsip wsink
243     & KsP KsN KsFe KsSi palat1 palat2 mQY(-4) chl2c Topt nsrc
244     & np'
245     do np = 1, npmax
246     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
247     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
248     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
249     & wsink(np),
250     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
251     & ,KsatSi(np),
252     & palat(np,1),palat(np,2),
253     & mQyield(np)*1e4,chl2cmax(np),
254     & phytoTempOptimum(np),nsource(np),np
255     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
256     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
257     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
258     & wsink(np),
259     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
260     & ,KsatSi(np),
261     & palat(np,1),palat(np,2),
262     & mQyield(np)*1e4,chl2cmax(np),
263     & phytoTempOptimum(np),nsource(np),np
264     end do
265     #endif
266    
267     close(IniUnit2)
268     close(IniUnit1)
269     110 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2e11.2,2f9.4,f6.1,2i5)
270     111 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2f6.1,2f9.4,f6.1,2i5)
271    
272     CALL LEF_ZERO( fice,myThid )
273     CALL LEF_ZERO( inputFe,myThid )
274     CALL LEF_ZERO( sur_par,myThid )
275     #ifdef NUT_SUPPLY
276     DO bj = myByLo(myThid), myByHi(myThid)
277     DO bi = myBxLo(myThid), myBxHi(myThid)
278     DO j=1-Oly,sNy+Oly
279     DO i=1-Olx,sNx+Olx
280     DO k=1,nR
281     nut_wvel(i,j,k,bi,bj) = 0. _d 0
282     ENDDO
283     ENDDO
284     ENDDO
285     ENDDO
286     ENDDO
287     #endif
288    
289     #ifdef ALLOW_PAR_DAY
290     DO iPAR=1,2
291     DO bj=myByLo(myThid), myByHi(myThid)
292     DO bi=myBxLo(myThid), myBxHi(myThid)
293     DO k=1,nR
294     DO j=1-Oly,sNy+Oly
295     DO i=1-Olx,sNx+Olx
296     PARday(i,j,k,bi,bj,iPAR) = 0. _d 0
297     ENDDO
298     ENDDO
299     ENDDO
300     ENDDO
301     ENDDO
302     ENDDO
303     IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
304     & .AND. pickupSuff .EQ. ' ') ) THEN
305     COJ should probably initialize from a file when nIter0 .EQ. 0
306     CALL DARWIN_READ_PICKUP( nIter0, myThid )
307     ENDIF
308     #endif
309     c
310     #ifdef ALLOW_TIMEAVE
311     c set arrays to zero if first timestep
312     DO bj = myByLo(myThid), myByHi(myThid)
313     DO bi = myBxLo(myThid), myBxHi(myThid)
314     CALL TIMEAVE_RESET(PARave, Nr, bi, bj, myThid)
315     CALL TIMEAVE_RESET(PPave, Nr, bi, bj, myThid)
316     CALL TIMEAVE_RESET(Chlave, Nr, bi, bj, myThid)
317     CALL TIMEAVE_RESET(Nfixave, Nr, bi, bj, myThid)
318     CALL TIMEAVE_RESET(Denitave, Nr, bi, bj, myThid)
319     #ifdef DAR_DIAG_ACDOM
320     CALL TIMEAVE_RESET(aCDOMave, Nr, bi, bj, myThid)
321     #endif
322     #ifdef DAR_DIAG_IRR
323     do i=1,tlam
324     CALL TIMEAVE_RESET(Edave(1-OLx,1-OLy,1,1,1,i),
325     & Nr,bi,bj,myThid)
326     CALL TIMEAVE_RESET(Esave(1-OLx,1-OLy,1,1,1,i),
327     & Nr,bi,bj,myThid)
328     CALL TIMEAVE_RESET(Euave(1-OLx,1-OLy,1,1,1,i),
329     & Nr,bi,bj,myThid)
330     CALL TIMEAVE_RESET(Eutave(1-OLx,1-OLy,1,1,1,i),
331     & Nr,bi,bj,myThid)
332     enddo
333     #endif
334     #ifdef DAR_DIAG_ABSORP
335     do i=1,tlam
336     CALL TIMEAVE_RESET(aave(1-OLx,1-OLy,1,1,1,i),
337     & Nr,bi,bj,myThid)
338     enddo
339     #endif
340     #ifdef DAR_DIAG_SCATTER
341     do i=1,tlam
342     CALL TIMEAVE_RESET(btave(1-OLx,1-OLy,1,1,1,i),
343     & Nr,bi,bj,myThid)
344     CALL TIMEAVE_RESET(bbave(1-OLx,1-OLy,1,1,1,i),
345     & Nr,bi,bj,myThid)
346     enddo
347     #endif
348     #ifdef DAR_DIAG_PART_SCATTER
349     do i=1,tlam
350     CALL TIMEAVE_RESET(apartave(1-OLx,1-OLy,1,1,1,i),
351     & Nr,bi,bj,myThid)
352     CALL TIMEAVE_RESET(btpartave(1-OLx,1-OLy,1,1,1,i),
353     & Nr,bi,bj,myThid)
354     CALL TIMEAVE_RESET(bbpartave(1-OLx,1-OLy,1,1,1,i),
355     & Nr,bi,bj,myThid)
356     enddo
357     #endif
358     c ANNA_TAVE
359     #ifdef WAVES_DIAG_PCHL
360     do np=1,npmax
361     CALL TIMEAVE_RESET(Pchlave(1-OLx,1-OLy,1,1,1,np),
362     & Nr,bi,bj,myThid)
363     enddo
364     #endif
365     c ANNA end TAVE
366     #ifdef DAR_DIAG_RSTAR
367     do np=1,npmax
368     CALL TIMEAVE_RESET(Rstarave(1-OLx,1-OLy,1,1,1,np),
369     & Nr,bi,bj,myThid)
370     CALL TIMEAVE_RESET(RNstarave(1-OLx,1-OLy,1,1,1,np),
371     & Nr,bi,bj,myThid)
372     enddo
373     #endif
374     #ifdef DAR_DIAG_DIVER
375     CALL TIMEAVE_RESET(Diver1ave, Nr, bi, bj, myThid)
376     CALL TIMEAVE_RESET(Diver2ave, Nr, bi, bj, myThid)
377     CALL TIMEAVE_RESET(Diver3ave, Nr, bi, bj, myThid)
378     CALL TIMEAVE_RESET(Diver4ave, Nr, bi, bj, myThid)
379     #endif
380     #ifdef ALLOW_DIAZ
381     #ifdef DAR_DIAG_NFIXP
382     do np=1,npmax
383     CALL TIMEAVE_RESET(NfixPave(1-OLx,1-OLy,1,1,1,np),
384     & Nr,bi,bj,myThid)
385     enddo
386     #endif
387     #endif
388     c CALL TIMEAVE_RESET(SURave, 1, bi, bj, myThid)
389     WRITE(msgbuf,'(A)')
390     & 'QQ start timeave'
391     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
392     & SQUEEZE_RIGHT , mythid)
393    
394     do k=1,Nr
395     DAR_TimeAve(bi,bj,k)=0. _d 0
396     enddo
397     ENDDO
398     ENDDO
399     #endif /* ALLOW_TIMEAVE */
400    
401     #ifdef CHECK_CONS
402     coj find unused units for darwin_cons output
403     CALL MDSFINDUNIT( DAR_cons_unit1, mythid )
404     open(DAR_cons_unit1,file='darwin_cons_P.txt',status='unknown')
405     CALL MDSFINDUNIT( DAR_cons_unit2, mythid )
406     open(DAR_cons_unit2,file='darwin_cons_N.txt',status='unknown')
407     CALL MDSFINDUNIT( DAR_cons_unit3, mythid )
408     open(DAR_cons_unit3,file='darwin_cons_Fe.txt',status='unknown')
409     CALL MDSFINDUNIT( DAR_cons_unit4, mythid )
410     open(DAR_cons_unit4,file='darwin_cons_Si.txt',status='unknown')
411     #ifdef ALLOW_CARBON
412     CALL MDSFINDUNIT( DAR_cons_unit5, mythid )
413     open(DAR_cons_unit5,file='darwin_cons_C.txt',status='unknown')
414     CALL MDSFINDUNIT( DAR_cons_unit6, mythid )
415     open(DAR_cons_unit6,file='darwin_cons_A.txt',status='unknown')
416     CALL MDSFINDUNIT( DAR_cons_unit7, mythid )
417     open(DAR_cons_unit7,file='darwin_cons_O.txt',status='unknown')
418     #endif
419     #endif
420     _END_MASTER(myThid)
421    
422     c test....................
423     c write(6,*)'finishing darwin_init_vari '
424     c test....................
425     WRITE(msgBuf,'(A)')
426     & '// ======================================================='
427     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
428     & SQUEEZE_RIGHT, myThid )
429     WRITE(msgBuf,'(A)') '// Darwin init variables >>> END <<<'
430     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
431     & SQUEEZE_RIGHT, myThid )
432     WRITE(msgBuf,'(A)')
433     & '// ======================================================='
434     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
435     & SQUEEZE_RIGHT, myThid )
436    
437    
438     RETURN
439     END
440     #endif /*MONOD*/
441     #endif /*ALLOW_PTRACERS*/
442     c ==========================================================
443    

  ViewVC Help
Powered by ViewVC 1.1.22