/[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.4 - (hide annotations) (download)
Wed Dec 7 20:04:00 2011 UTC (13 years, 7 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63p_20120707, ctrb_darwin2_ckpt63n_20120604, ctrb_darwin2_ckpt63k_20120317
Changes since 1.3: +13 -3 lines
fix some multi-threading bits

1 jahn 1.4 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_init_vari.F,v 1.3 2011/09/22 16:01:30 jahn 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 jahn 1.4 _BEGIN_MASTER( myThid )
92    
93 jahn 1.1 CALL DARWIN_RANDOM_INIT(darwin_seed, myThid)
94    
95     c initialize total number of functional groups tried
96     ngroups = 0
97     do np = 1, npmax
98     #ifdef ALLOW_MUTANTS
99     call MONOD_GENERATE_MUTANTS(MyThid, np)
100     #else
101     call MONOD_GENERATE_PHYTO(MyThid, np)
102     #endif
103     end do
104    
105 jahn 1.4 _END_MASTER( myThid )
106    
107 jahn 1.1 c reduce amount of diaz
108     #ifdef ALLOW_DIAZ
109     do np = 1, npmax
110     if (diazotroph(np) .eq. 1. _d 0) then
111     DO bj = myByLo(myThid), myByHi(myThid)
112     DO bi = myBxLo(myThid), myBxHi(myThid)
113     DO j=1-Oly,sNy+Oly
114     DO i=1-Olx,sNx+Olx
115     DO k=1,nR
116     Ptracer(i,j,k,bi,bj,iPhy+np-1) =
117     & Ptracer(i,j,k,bi,bj,iPhy+np-1)/10. _d 0
118     ENDDO
119     ENDDO
120     ENDDO
121     ENDDO
122     ENDDO
123     endif
124     enddo
125     #endif
126    
127 jahn 1.4 _BEGIN_MASTER( myThid )
128    
129 jahn 1.1 c initialize zooplankton
130     call MONOD_GENERATE_ZOO(MyThid)
131    
132     c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133    
134     c ANNA call WAVEBANDS_INIT_VARI to assign variable parameters
135     #ifdef WAVEBANDS
136     call WAVEBANDS_INIT_VARI(MyThid)
137     #endif
138    
139     c ANNA get alphachl from mQyield / aphy_chl
140     c ANNA must do this after params are assigned, but before written out
141     c ANNA use aphy_chl_ps for growth. To turn off, simply set same coefs as aphy_chl in input files.
142     #ifdef GEIDER
143     #ifndef WAVEBANDS
144     do np = 1,npmax
145     alphachl(np) = mQyield(np) * aphy_chl_ave
146     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
147     c no temp/nutrient limitation
148     c chl2cmin(np)=chl2cmax(np)/
149     c & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
150     c & (2*pcmax(np)))
151     chl2cmin(np)=0. _d 0
152     enddo
153     #else
154     do np = 1,npmax
155     do nl = 1,tlam
156     alphachl_nl(np,nl) = mQyield(np) * aphy_chl_ps(np,nl)
157     end do
158     c find mean
159     cu_area = 0.d0
160     do nl = 1,tlam
161     cu_area = cu_area + wb_width(nl)*alphachl_nl(np,nl)
162     end do
163     alpha_mean(np) = cu_area / wb_totalWidth
164    
165     chl2cmin(np)=chl2cmax(np)/
166     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
167     & (2*pcmax(np)))
168     end do
169     #endif
170     #ifdef DYNAMIC_CHL
171     c check Chl fields are reasonable
172 stephd 1.2 #ifndef WAVEBANDS
173     do np = 1,npmax
174     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
175     c no temp/nutrient limitation
176     chl2cmin(np)=chl2cmax(np)/
177     & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
178     & (2*pcmax(np)))
179     chl2cmin(np)=0. _d 0
180     enddo
181     #else
182 jahn 1.1 do np=1,npmax
183     chl2cmin(np)=chl2cmax(np)/
184     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
185     & (2*pcmax(np)))
186     enddo
187 stephd 1.2 #endif
188 jahn 1.1 #endif
189     #endif
190 jahn 1.4 _END_MASTER( myThid )
191    
192     #ifdef DYNAMIC_CHL
193     C this initializes fields...
194     call MONOD_CHECK_CHL(myThid)
195     #endif
196 jahn 1.1
197 jahn 1.3 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
198 jahn 1.1 c write out initial phyto characteristics
199     #ifndef GEIDER
200     CALL MDSFINDUNIT( IniUnit1, mythid )
201     open(IniUnit1,file='plankton-ini-char.dat',status='unknown')
202     CALL MDSFINDUNIT( IniUnit2, mythid )
203     open(IniUnit2,file='plankton_ini_char_nohead.dat',
204     & status='unknown')
205     #ifdef OLD_GRAZE
206     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip
207     & wsink KsP KsN KsFe KsSi g1 g2 Kpar Kinh
208     & Topt nsrc np'
209     do np = 1, npmax
210     write(IniUnit1,110)diacoc(np),diazotroph(np),physize(np),
211     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
212     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
213     & wsink(np),
214     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
215     & ,KsatSi(np),
216     & graze(np,1),graze(np,2),
217     & KsatPAR(np),Kinhib(np),
218     & phytoTempOptimum(np),nsource(np),np
219     write(IniUnit2,110)diacoc(np),diazotroph(np),physize(np),
220     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
221     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
222     & wsink(np),
223     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
224     & ,KsatSi(np),
225     & graze(np,1),graze(np,2),
226     & KsatPAR(np),Kinhib(np),
227     & phytoTempOptimum(np),nsource(np),np
228     end do
229     #else
230     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip wsink
231     & KsP KsN KsFe KsSi palat1 palat2 Kpar Kinh Topt nsrc
232     & np'
233     do np = 1, npmax
234     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
235     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
236     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
237     & wsink(np),
238     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
239     & ,KsatSi(np),
240     & palat(np,1),palat(np,2),
241     & KsatPAR(np),Kinhib(np),
242     & phytoTempOptimum(np),nsource(np),np
243     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
244     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
245     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
246     & wsink(np),
247     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
248     & ,KsatSi(np),
249     & palat(np,1),palat(np,2),
250     & KsatPAR(np),Kinhib(np),
251     & phytoTempOptimum(np),nsource(np),np
252     end do
253     #endif
254     #endif
255    
256     #ifdef GEIDER
257     c ANNA outputs mQyield as 10^(4) mmol C (uEin)-1
258     CALL MDSFINDUNIT( IniUnit1, mythid )
259     open(IniUnit1,file='gplankton-ini-char.dat',status='unknown')
260     CALL MDSFINDUNIT( IniUnit2, mythid )
261     open(IniUnit2,file='gplankton_ini_char_nohead.dat',
262     & status='unknown')
263     write(IniUnit1,*)'dico diaz size pcmax mort Rnp Rfep Rsip wsink
264     & KsP KsN KsFe KsSi palat1 palat2 mQY(-4) chl2c Topt nsrc
265     & np'
266     do np = 1, npmax
267     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
268     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
269     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
270     & wsink(np),
271     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
272     & ,KsatSi(np),
273     & palat(np,1),palat(np,2),
274     & mQyield(np)*1e4,chl2cmax(np),
275     & phytoTempOptimum(np),nsource(np),np
276     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
277     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
278     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
279     & wsink(np),
280     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
281     & ,KsatSi(np),
282     & palat(np,1),palat(np,2),
283     & mQyield(np)*1e4,chl2cmax(np),
284     & phytoTempOptimum(np),nsource(np),np
285     end do
286     #endif
287    
288     close(IniUnit2)
289     close(IniUnit1)
290     110 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2e11.2,2f9.4,f6.1,2i5)
291     111 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2f6.1,2f9.4,f6.1,2i5)
292 jahn 1.3 c myProcId and myThid
293     ENDIF
294 jahn 1.1
295     CALL LEF_ZERO( fice,myThid )
296     CALL LEF_ZERO( inputFe,myThid )
297     CALL LEF_ZERO( sur_par,myThid )
298     #ifdef NUT_SUPPLY
299     DO bj = myByLo(myThid), myByHi(myThid)
300     DO bi = myBxLo(myThid), myBxHi(myThid)
301     DO j=1-Oly,sNy+Oly
302     DO i=1-Olx,sNx+Olx
303     DO k=1,nR
304     nut_wvel(i,j,k,bi,bj) = 0. _d 0
305     ENDDO
306     ENDDO
307     ENDDO
308     ENDDO
309     ENDDO
310     #endif
311    
312     #ifdef ALLOW_PAR_DAY
313     DO iPAR=1,2
314     DO bj=myByLo(myThid), myByHi(myThid)
315     DO bi=myBxLo(myThid), myBxHi(myThid)
316     DO k=1,nR
317     DO j=1-Oly,sNy+Oly
318     DO i=1-Olx,sNx+Olx
319     PARday(i,j,k,bi,bj,iPAR) = 0. _d 0
320     ENDDO
321     ENDDO
322     ENDDO
323     ENDDO
324     ENDDO
325     ENDDO
326     IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
327     & .AND. pickupSuff .EQ. ' ') ) THEN
328     COJ should probably initialize from a file when nIter0 .EQ. 0
329     CALL DARWIN_READ_PICKUP( nIter0, myThid )
330     ENDIF
331     #endif
332     c
333     #ifdef ALLOW_TIMEAVE
334     c set arrays to zero if first timestep
335     DO bj = myByLo(myThid), myByHi(myThid)
336     DO bi = myBxLo(myThid), myBxHi(myThid)
337     CALL TIMEAVE_RESET(PARave, Nr, bi, bj, myThid)
338     CALL TIMEAVE_RESET(PPave, Nr, bi, bj, myThid)
339     CALL TIMEAVE_RESET(Chlave, Nr, bi, bj, myThid)
340     CALL TIMEAVE_RESET(Nfixave, Nr, bi, bj, myThid)
341     CALL TIMEAVE_RESET(Denitave, Nr, bi, bj, myThid)
342     #ifdef DAR_DIAG_ACDOM
343     CALL TIMEAVE_RESET(aCDOMave, Nr, bi, bj, myThid)
344     #endif
345     #ifdef DAR_DIAG_IRR
346     do i=1,tlam
347     CALL TIMEAVE_RESET(Edave(1-OLx,1-OLy,1,1,1,i),
348     & Nr,bi,bj,myThid)
349     CALL TIMEAVE_RESET(Esave(1-OLx,1-OLy,1,1,1,i),
350     & Nr,bi,bj,myThid)
351     CALL TIMEAVE_RESET(Euave(1-OLx,1-OLy,1,1,1,i),
352     & Nr,bi,bj,myThid)
353     CALL TIMEAVE_RESET(Eutave(1-OLx,1-OLy,1,1,1,i),
354     & Nr,bi,bj,myThid)
355     enddo
356     #endif
357     #ifdef DAR_DIAG_ABSORP
358     do i=1,tlam
359     CALL TIMEAVE_RESET(aave(1-OLx,1-OLy,1,1,1,i),
360     & Nr,bi,bj,myThid)
361     enddo
362     #endif
363     #ifdef DAR_DIAG_SCATTER
364     do i=1,tlam
365     CALL TIMEAVE_RESET(btave(1-OLx,1-OLy,1,1,1,i),
366     & Nr,bi,bj,myThid)
367     CALL TIMEAVE_RESET(bbave(1-OLx,1-OLy,1,1,1,i),
368     & Nr,bi,bj,myThid)
369     enddo
370     #endif
371     #ifdef DAR_DIAG_PART_SCATTER
372     do i=1,tlam
373     CALL TIMEAVE_RESET(apartave(1-OLx,1-OLy,1,1,1,i),
374     & Nr,bi,bj,myThid)
375     CALL TIMEAVE_RESET(btpartave(1-OLx,1-OLy,1,1,1,i),
376     & Nr,bi,bj,myThid)
377     CALL TIMEAVE_RESET(bbpartave(1-OLx,1-OLy,1,1,1,i),
378     & Nr,bi,bj,myThid)
379     enddo
380     #endif
381     c ANNA_TAVE
382     #ifdef WAVES_DIAG_PCHL
383     do np=1,npmax
384     CALL TIMEAVE_RESET(Pchlave(1-OLx,1-OLy,1,1,1,np),
385     & Nr,bi,bj,myThid)
386     enddo
387     #endif
388     c ANNA end TAVE
389     #ifdef DAR_DIAG_RSTAR
390     do np=1,npmax
391     CALL TIMEAVE_RESET(Rstarave(1-OLx,1-OLy,1,1,1,np),
392     & Nr,bi,bj,myThid)
393     CALL TIMEAVE_RESET(RNstarave(1-OLx,1-OLy,1,1,1,np),
394     & Nr,bi,bj,myThid)
395     enddo
396     #endif
397     #ifdef DAR_DIAG_DIVER
398     CALL TIMEAVE_RESET(Diver1ave, Nr, bi, bj, myThid)
399     CALL TIMEAVE_RESET(Diver2ave, Nr, bi, bj, myThid)
400     CALL TIMEAVE_RESET(Diver3ave, Nr, bi, bj, myThid)
401     CALL TIMEAVE_RESET(Diver4ave, Nr, bi, bj, myThid)
402     #endif
403     #ifdef ALLOW_DIAZ
404     #ifdef DAR_DIAG_NFIXP
405     do np=1,npmax
406     CALL TIMEAVE_RESET(NfixPave(1-OLx,1-OLy,1,1,1,np),
407     & Nr,bi,bj,myThid)
408     enddo
409     #endif
410     #endif
411     c CALL TIMEAVE_RESET(SURave, 1, bi, bj, myThid)
412     WRITE(msgbuf,'(A)')
413     & 'QQ start timeave'
414     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
415     & SQUEEZE_RIGHT , mythid)
416    
417     do k=1,Nr
418     DAR_TimeAve(bi,bj,k)=0. _d 0
419     enddo
420     ENDDO
421     ENDDO
422     #endif /* ALLOW_TIMEAVE */
423    
424     #ifdef CHECK_CONS
425 jahn 1.3 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
426 jahn 1.1 coj find unused units for darwin_cons output
427     CALL MDSFINDUNIT( DAR_cons_unit1, mythid )
428     open(DAR_cons_unit1,file='darwin_cons_P.txt',status='unknown')
429     CALL MDSFINDUNIT( DAR_cons_unit2, mythid )
430     open(DAR_cons_unit2,file='darwin_cons_N.txt',status='unknown')
431     CALL MDSFINDUNIT( DAR_cons_unit3, mythid )
432     open(DAR_cons_unit3,file='darwin_cons_Fe.txt',status='unknown')
433     CALL MDSFINDUNIT( DAR_cons_unit4, mythid )
434     open(DAR_cons_unit4,file='darwin_cons_Si.txt',status='unknown')
435     #ifdef ALLOW_CARBON
436     CALL MDSFINDUNIT( DAR_cons_unit5, mythid )
437     open(DAR_cons_unit5,file='darwin_cons_C.txt',status='unknown')
438     CALL MDSFINDUNIT( DAR_cons_unit6, mythid )
439     open(DAR_cons_unit6,file='darwin_cons_A.txt',status='unknown')
440     CALL MDSFINDUNIT( DAR_cons_unit7, mythid )
441     open(DAR_cons_unit7,file='darwin_cons_O.txt',status='unknown')
442     #endif
443 jahn 1.3 c myProcId and myThid
444     ENDIF
445 jahn 1.1 #endif
446    
447     c test....................
448     c write(6,*)'finishing darwin_init_vari '
449     c test....................
450     WRITE(msgBuf,'(A)')
451     & '// ======================================================='
452     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
453     & SQUEEZE_RIGHT, myThid )
454     WRITE(msgBuf,'(A)') '// Darwin init variables >>> END <<<'
455     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
456     & SQUEEZE_RIGHT, myThid )
457     WRITE(msgBuf,'(A)')
458     & '// ======================================================='
459     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
460     & SQUEEZE_RIGHT, myThid )
461    
462    
463     RETURN
464     END
465     #endif /*MONOD*/
466     #endif /*ALLOW_PTRACERS*/
467     c ==========================================================
468    

  ViewVC Help
Powered by ViewVC 1.1.22