/[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.3 - (hide annotations) (download)
Thu Sep 22 16:01:30 2011 UTC (13 years, 9 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt63c_20111011, ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63d_20111107
Changes since 1.2: +7 -3 lines
only master cpu/thread writes

1 jahn 1.3 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_init_vari.F,v 1.2 2011/05/02 17:55:00 stephd Exp $
2 stephd 1.2 C $Name: $
3 jahn 1.1
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 stephd 1.2 #ifndef WAVEBANDS
167     do np = 1,npmax
168     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
169     c no temp/nutrient limitation
170     chl2cmin(np)=chl2cmax(np)/
171     & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
172     & (2*pcmax(np)))
173     chl2cmin(np)=0. _d 0
174     enddo
175     #else
176 jahn 1.1 do np=1,npmax
177     chl2cmin(np)=chl2cmax(np)/
178     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
179     & (2*pcmax(np)))
180     enddo
181 stephd 1.2 #endif
182 jahn 1.1 call MONOD_CHECK_CHL(myThid)
183     #endif
184     #endif
185     c ANNA endif
186    
187 jahn 1.3 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
188 jahn 1.1 c write out initial phyto characteristics
189     #ifndef GEIDER
190     CALL MDSFINDUNIT( IniUnit1, mythid )
191     open(IniUnit1,file='plankton-ini-char.dat',status='unknown')
192     CALL MDSFINDUNIT( IniUnit2, mythid )
193     open(IniUnit2,file='plankton_ini_char_nohead.dat',
194     & status='unknown')
195     #ifdef OLD_GRAZE
196     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip
197     & wsink KsP KsN KsFe KsSi g1 g2 Kpar Kinh
198     & Topt nsrc np'
199     do np = 1, npmax
200     write(IniUnit1,110)diacoc(np),diazotroph(np),physize(np),
201     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
202     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
203     & wsink(np),
204     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
205     & ,KsatSi(np),
206     & graze(np,1),graze(np,2),
207     & KsatPAR(np),Kinhib(np),
208     & phytoTempOptimum(np),nsource(np),np
209     write(IniUnit2,110)diacoc(np),diazotroph(np),physize(np),
210     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
211     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
212     & wsink(np),
213     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
214     & ,KsatSi(np),
215     & graze(np,1),graze(np,2),
216     & KsatPAR(np),Kinhib(np),
217     & phytoTempOptimum(np),nsource(np),np
218     end do
219     #else
220     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip wsink
221     & KsP KsN KsFe KsSi palat1 palat2 Kpar Kinh Topt nsrc
222     & np'
223     do np = 1, npmax
224     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
225     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
226     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
227     & wsink(np),
228     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
229     & ,KsatSi(np),
230     & palat(np,1),palat(np,2),
231     & KsatPAR(np),Kinhib(np),
232     & phytoTempOptimum(np),nsource(np),np
233     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
234     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
235     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
236     & wsink(np),
237     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
238     & ,KsatSi(np),
239     & palat(np,1),palat(np,2),
240     & KsatPAR(np),Kinhib(np),
241     & phytoTempOptimum(np),nsource(np),np
242     end do
243     #endif
244     #endif
245    
246     #ifdef GEIDER
247     c ANNA outputs mQyield as 10^(4) mmol C (uEin)-1
248     CALL MDSFINDUNIT( IniUnit1, mythid )
249     open(IniUnit1,file='gplankton-ini-char.dat',status='unknown')
250     CALL MDSFINDUNIT( IniUnit2, mythid )
251     open(IniUnit2,file='gplankton_ini_char_nohead.dat',
252     & status='unknown')
253     write(IniUnit1,*)'dico diaz size pcmax mort Rnp Rfep Rsip wsink
254     & KsP KsN KsFe KsSi palat1 palat2 mQY(-4) chl2c Topt nsrc
255     & np'
256     do np = 1, npmax
257     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
258     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
259     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
260     & wsink(np),
261     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
262     & ,KsatSi(np),
263     & palat(np,1),palat(np,2),
264     & mQyield(np)*1e4,chl2cmax(np),
265     & phytoTempOptimum(np),nsource(np),np
266     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
267     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
268     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
269     & wsink(np),
270     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
271     & ,KsatSi(np),
272     & palat(np,1),palat(np,2),
273     & mQyield(np)*1e4,chl2cmax(np),
274     & phytoTempOptimum(np),nsource(np),np
275     end do
276     #endif
277    
278     close(IniUnit2)
279     close(IniUnit1)
280     110 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2e11.2,2f9.4,f6.1,2i5)
281     111 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2f6.1,2f9.4,f6.1,2i5)
282 jahn 1.3 c myProcId and myThid
283     ENDIF
284 jahn 1.1
285     CALL LEF_ZERO( fice,myThid )
286     CALL LEF_ZERO( inputFe,myThid )
287     CALL LEF_ZERO( sur_par,myThid )
288     #ifdef NUT_SUPPLY
289     DO bj = myByLo(myThid), myByHi(myThid)
290     DO bi = myBxLo(myThid), myBxHi(myThid)
291     DO j=1-Oly,sNy+Oly
292     DO i=1-Olx,sNx+Olx
293     DO k=1,nR
294     nut_wvel(i,j,k,bi,bj) = 0. _d 0
295     ENDDO
296     ENDDO
297     ENDDO
298     ENDDO
299     ENDDO
300     #endif
301    
302     #ifdef ALLOW_PAR_DAY
303     DO iPAR=1,2
304     DO bj=myByLo(myThid), myByHi(myThid)
305     DO bi=myBxLo(myThid), myBxHi(myThid)
306     DO k=1,nR
307     DO j=1-Oly,sNy+Oly
308     DO i=1-Olx,sNx+Olx
309     PARday(i,j,k,bi,bj,iPAR) = 0. _d 0
310     ENDDO
311     ENDDO
312     ENDDO
313     ENDDO
314     ENDDO
315     ENDDO
316     IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
317     & .AND. pickupSuff .EQ. ' ') ) THEN
318     COJ should probably initialize from a file when nIter0 .EQ. 0
319     CALL DARWIN_READ_PICKUP( nIter0, myThid )
320     ENDIF
321     #endif
322     c
323     #ifdef ALLOW_TIMEAVE
324     c set arrays to zero if first timestep
325     DO bj = myByLo(myThid), myByHi(myThid)
326     DO bi = myBxLo(myThid), myBxHi(myThid)
327     CALL TIMEAVE_RESET(PARave, Nr, bi, bj, myThid)
328     CALL TIMEAVE_RESET(PPave, Nr, bi, bj, myThid)
329     CALL TIMEAVE_RESET(Chlave, Nr, bi, bj, myThid)
330     CALL TIMEAVE_RESET(Nfixave, Nr, bi, bj, myThid)
331     CALL TIMEAVE_RESET(Denitave, Nr, bi, bj, myThid)
332     #ifdef DAR_DIAG_ACDOM
333     CALL TIMEAVE_RESET(aCDOMave, Nr, bi, bj, myThid)
334     #endif
335     #ifdef DAR_DIAG_IRR
336     do i=1,tlam
337     CALL TIMEAVE_RESET(Edave(1-OLx,1-OLy,1,1,1,i),
338     & Nr,bi,bj,myThid)
339     CALL TIMEAVE_RESET(Esave(1-OLx,1-OLy,1,1,1,i),
340     & Nr,bi,bj,myThid)
341     CALL TIMEAVE_RESET(Euave(1-OLx,1-OLy,1,1,1,i),
342     & Nr,bi,bj,myThid)
343     CALL TIMEAVE_RESET(Eutave(1-OLx,1-OLy,1,1,1,i),
344     & Nr,bi,bj,myThid)
345     enddo
346     #endif
347     #ifdef DAR_DIAG_ABSORP
348     do i=1,tlam
349     CALL TIMEAVE_RESET(aave(1-OLx,1-OLy,1,1,1,i),
350     & Nr,bi,bj,myThid)
351     enddo
352     #endif
353     #ifdef DAR_DIAG_SCATTER
354     do i=1,tlam
355     CALL TIMEAVE_RESET(btave(1-OLx,1-OLy,1,1,1,i),
356     & Nr,bi,bj,myThid)
357     CALL TIMEAVE_RESET(bbave(1-OLx,1-OLy,1,1,1,i),
358     & Nr,bi,bj,myThid)
359     enddo
360     #endif
361     #ifdef DAR_DIAG_PART_SCATTER
362     do i=1,tlam
363     CALL TIMEAVE_RESET(apartave(1-OLx,1-OLy,1,1,1,i),
364     & Nr,bi,bj,myThid)
365     CALL TIMEAVE_RESET(btpartave(1-OLx,1-OLy,1,1,1,i),
366     & Nr,bi,bj,myThid)
367     CALL TIMEAVE_RESET(bbpartave(1-OLx,1-OLy,1,1,1,i),
368     & Nr,bi,bj,myThid)
369     enddo
370     #endif
371     c ANNA_TAVE
372     #ifdef WAVES_DIAG_PCHL
373     do np=1,npmax
374     CALL TIMEAVE_RESET(Pchlave(1-OLx,1-OLy,1,1,1,np),
375     & Nr,bi,bj,myThid)
376     enddo
377     #endif
378     c ANNA end TAVE
379     #ifdef DAR_DIAG_RSTAR
380     do np=1,npmax
381     CALL TIMEAVE_RESET(Rstarave(1-OLx,1-OLy,1,1,1,np),
382     & Nr,bi,bj,myThid)
383     CALL TIMEAVE_RESET(RNstarave(1-OLx,1-OLy,1,1,1,np),
384     & Nr,bi,bj,myThid)
385     enddo
386     #endif
387     #ifdef DAR_DIAG_DIVER
388     CALL TIMEAVE_RESET(Diver1ave, Nr, bi, bj, myThid)
389     CALL TIMEAVE_RESET(Diver2ave, Nr, bi, bj, myThid)
390     CALL TIMEAVE_RESET(Diver3ave, Nr, bi, bj, myThid)
391     CALL TIMEAVE_RESET(Diver4ave, Nr, bi, bj, myThid)
392     #endif
393     #ifdef ALLOW_DIAZ
394     #ifdef DAR_DIAG_NFIXP
395     do np=1,npmax
396     CALL TIMEAVE_RESET(NfixPave(1-OLx,1-OLy,1,1,1,np),
397     & Nr,bi,bj,myThid)
398     enddo
399     #endif
400     #endif
401     c CALL TIMEAVE_RESET(SURave, 1, bi, bj, myThid)
402     WRITE(msgbuf,'(A)')
403     & 'QQ start timeave'
404     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
405     & SQUEEZE_RIGHT , mythid)
406    
407     do k=1,Nr
408     DAR_TimeAve(bi,bj,k)=0. _d 0
409     enddo
410     ENDDO
411     ENDDO
412     #endif /* ALLOW_TIMEAVE */
413    
414     #ifdef CHECK_CONS
415 jahn 1.3 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
416 jahn 1.1 coj find unused units for darwin_cons output
417     CALL MDSFINDUNIT( DAR_cons_unit1, mythid )
418     open(DAR_cons_unit1,file='darwin_cons_P.txt',status='unknown')
419     CALL MDSFINDUNIT( DAR_cons_unit2, mythid )
420     open(DAR_cons_unit2,file='darwin_cons_N.txt',status='unknown')
421     CALL MDSFINDUNIT( DAR_cons_unit3, mythid )
422     open(DAR_cons_unit3,file='darwin_cons_Fe.txt',status='unknown')
423     CALL MDSFINDUNIT( DAR_cons_unit4, mythid )
424     open(DAR_cons_unit4,file='darwin_cons_Si.txt',status='unknown')
425     #ifdef ALLOW_CARBON
426     CALL MDSFINDUNIT( DAR_cons_unit5, mythid )
427     open(DAR_cons_unit5,file='darwin_cons_C.txt',status='unknown')
428     CALL MDSFINDUNIT( DAR_cons_unit6, mythid )
429     open(DAR_cons_unit6,file='darwin_cons_A.txt',status='unknown')
430     CALL MDSFINDUNIT( DAR_cons_unit7, mythid )
431     open(DAR_cons_unit7,file='darwin_cons_O.txt',status='unknown')
432     #endif
433 jahn 1.3 c myProcId and myThid
434     ENDIF
435 jahn 1.1 #endif
436    
437     c test....................
438     c write(6,*)'finishing darwin_init_vari '
439     c test....................
440     WRITE(msgBuf,'(A)')
441     & '// ======================================================='
442     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
443     & SQUEEZE_RIGHT, myThid )
444     WRITE(msgBuf,'(A)') '// Darwin init variables >>> END <<<'
445     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
446     & SQUEEZE_RIGHT, myThid )
447     WRITE(msgBuf,'(A)')
448     & '// ======================================================='
449     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
450     & SQUEEZE_RIGHT, myThid )
451    
452    
453     RETURN
454     END
455     #endif /*MONOD*/
456     #endif /*ALLOW_PTRACERS*/
457     c ==========================================================
458    

  ViewVC Help
Powered by ViewVC 1.1.22