/[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.13 - (hide annotations) (download)
Wed Dec 4 21:21:33 2013 UTC (11 years, 7 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64r_20131210
Changes since 1.12: +80 -92 lines
clean up multi-threading directives

1 jahn 1.13 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_init_vari.F,v 1.12 2013/06/21 12:49:47 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 jahn 1.13 C !INPUT PARAMETERS: ===================================================
45     C myThid :: thread number
46 jahn 1.1 INTEGER myThid
47 jahn 1.13 CEOP
48 jahn 1.1
49 jahn 1.13 C !FUNCTIONS:
50 jahn 1.1 _RL DARWIN_RANDOM
51     EXTERNAL DARWIN_RANDOM
52    
53     C !LOCAL VARIABLES:
54     C msgBuf - Informational/error meesage buffer
55     CHARACTER*(MAX_LEN_MBUF) msgBuf
56     INTEGER IniUnit1, IniUnit2
57    
58 jahn 1.13 INTEGER bi, bj, k, i, j, iPAR
59 jahn 1.1
60 jahn 1.13 INTEGER np
61     INTEGER nz
62 jahn 1.1 c ANNA need nl for wavebands
63     #ifdef WAVEBANDS
64 jahn 1.13 INTEGER ilam
65     INTEGER nl
66     _RL cu_area
67 jahn 1.1 #endif
68    
69 jahn 1.13 C ----------------------------------------------------------------------
70     C Scalar bits first
71     _BEGIN_MASTER( myThid )
72 jahn 1.1
73     WRITE(msgBuf,'(A)')
74     & '// ======================================================='
75     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
76     & SQUEEZE_RIGHT, myThid )
77     WRITE(msgBuf,'(A)') '// Darwin init variables >>> START <<<'
78     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
79     & SQUEEZE_RIGHT, myThid )
80     WRITE(msgBuf,'(A)')
81     & '// ======================================================='
82     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
83     & SQUEEZE_RIGHT, myThid )
84    
85     c set up ecosystem coefficients
86     c
87     c seed randomization
88 jahn 1.4
89 jahn 1.1 CALL DARWIN_RANDOM_INIT(darwin_seed, myThid)
90    
91     c initialize total number of functional groups tried
92     ngroups = 0
93     do np = 1, npmax
94     #ifdef ALLOW_MUTANTS
95     call MONOD_GENERATE_MUTANTS(MyThid, np)
96     #else
97     call MONOD_GENERATE_PHYTO(MyThid, np)
98     #endif
99     end do
100    
101     c initialize zooplankton
102     call MONOD_GENERATE_ZOO(MyThid)
103    
104     c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105    
106     c ANNA call WAVEBANDS_INIT_VARI to assign variable parameters
107     #ifdef WAVEBANDS
108     call WAVEBANDS_INIT_VARI(MyThid)
109     #endif
110    
111     c ANNA get alphachl from mQyield / aphy_chl
112     c ANNA must do this after params are assigned, but before written out
113     c ANNA use aphy_chl_ps for growth. To turn off, simply set same coefs as aphy_chl in input files.
114     #ifdef GEIDER
115     #ifndef WAVEBANDS
116     do np = 1,npmax
117     alphachl(np) = mQyield(np) * aphy_chl_ave
118     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
119     c no temp/nutrient limitation
120     c chl2cmin(np)=chl2cmax(np)/
121     c & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
122     c & (2*pcmax(np)))
123     chl2cmin(np)=0. _d 0
124     enddo
125     #else
126     do np = 1,npmax
127     do nl = 1,tlam
128     alphachl_nl(np,nl) = mQyield(np) * aphy_chl_ps(np,nl)
129     end do
130     c find mean
131     cu_area = 0.d0
132     do nl = 1,tlam
133     cu_area = cu_area + wb_width(nl)*alphachl_nl(np,nl)
134     end do
135     alpha_mean(np) = cu_area / wb_totalWidth
136    
137     chl2cmin(np)=chl2cmax(np)/
138     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
139     & (2*pcmax(np)))
140     end do
141     #endif
142     #ifdef DYNAMIC_CHL
143     c check Chl fields are reasonable
144 stephd 1.2 #ifndef WAVEBANDS
145     do np = 1,npmax
146     c C:CHl minimum: chosen to be Chl:C at high light (2000uEin/m2/s) and
147     c no temp/nutrient limitation
148     chl2cmin(np)=chl2cmax(np)/
149     & (1+(chl2cmax(np)*alphachl(np)*2000. _d 0)/
150     & (2*pcmax(np)))
151     chl2cmin(np)=0. _d 0
152     enddo
153     #else
154 jahn 1.1 do np=1,npmax
155     chl2cmin(np)=chl2cmax(np)/
156     & (1+(chl2cmax(np)* alpha_mean(np) *2000. _d 0)/
157     & (2*pcmax(np)))
158     enddo
159 stephd 1.2 #endif
160 jahn 1.1 #endif
161     #endif
162    
163 jahn 1.3 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
164 jahn 1.1 c write out initial phyto characteristics
165     #ifndef GEIDER
166     CALL MDSFINDUNIT( IniUnit1, mythid )
167     open(IniUnit1,file='plankton-ini-char.dat',status='unknown')
168     CALL MDSFINDUNIT( IniUnit2, mythid )
169     open(IniUnit2,file='plankton_ini_char_nohead.dat',
170     & status='unknown')
171     #ifdef OLD_GRAZE
172     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip
173     & wsink KsP KsN KsFe KsSi g1 g2 Kpar Kinh
174     & Topt nsrc np'
175     do np = 1, npmax
176     write(IniUnit1,110)diacoc(np),diazotroph(np),physize(np),
177     & 1.0/(mu(np)*86400.), 1.0/(mortphy(np)*86400.),
178     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
179     & wsink(np),
180     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
181     & ,KsatSi(np),
182     & graze(np,1),graze(np,2),
183     & KsatPAR(np),Kinhib(np),
184     & phytoTempOptimum(np),nsource(np),np
185     write(IniUnit2,110)diacoc(np),diazotroph(np),physize(np),
186     & 1.0/(mu(np)*86400.),1.0/(mortphy(np)*86400.),
187     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
188     & wsink(np),
189     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
190     & ,KsatSi(np),
191     & graze(np,1),graze(np,2),
192     & KsatPAR(np),Kinhib(np),
193     & phytoTempOptimum(np),nsource(np),np
194     end do
195     #else
196     write(IniUnit1,*)'dico diaz size mu mort Rnp Rfep Rsip wsink
197     & KsP KsN KsFe KsSi palat1 palat2 Kpar Kinh Topt nsrc
198     & np'
199     do np = 1, npmax
200     write(IniUnit1,111)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     & palat(np,1),palat(np,2),
207     & KsatPAR(np),Kinhib(np),
208     & phytoTempOptimum(np),nsource(np),np
209     write(IniUnit2,111)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     & palat(np,1),palat(np,2),
216     & KsatPAR(np),Kinhib(np),
217     & phytoTempOptimum(np),nsource(np),np
218     end do
219     #endif
220     #endif
221    
222     #ifdef GEIDER
223     c ANNA outputs mQyield as 10^(4) mmol C (uEin)-1
224     CALL MDSFINDUNIT( IniUnit1, mythid )
225     open(IniUnit1,file='gplankton-ini-char.dat',status='unknown')
226     CALL MDSFINDUNIT( IniUnit2, mythid )
227     open(IniUnit2,file='gplankton_ini_char_nohead.dat',
228     & status='unknown')
229     write(IniUnit1,*)'dico diaz size pcmax mort Rnp Rfep Rsip wsink
230     & KsP KsN KsFe KsSi palat1 palat2 mQY(-4) chl2c Topt nsrc
231     & np'
232     do np = 1, npmax
233     write(IniUnit1,111)diacoc(np),diazotroph(np),physize(np),
234     & 1.0/(pcmax(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     & mQyield(np)*1e4,chl2cmax(np),
241     & phytoTempOptimum(np),nsource(np),np
242     write(IniUnit2,111)diacoc(np),diazotroph(np),physize(np),
243     & 1.0/(pcmax(np)*86400.), 1.0/(mortphy(np)*86400.),
244     & R_NP(np),R_FeP(np)*1000.,R_SiP(np),
245     & wsink(np),
246     & KsatPO4(np),KsatNO3(np),KsatFeT(np)*1000.
247     & ,KsatSi(np),
248     & palat(np,1),palat(np,2),
249     & mQyield(np)*1e4,chl2cmax(np),
250     & phytoTempOptimum(np),nsource(np),np
251     end do
252     #endif
253    
254     close(IniUnit2)
255     close(IniUnit1)
256     110 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2e11.2,2f9.4,f6.1,2i5)
257     111 format(3f4.0,f6.2,4f4.0,f5.1,4f7.3,2f6.1,2f9.4,f6.1,2i5)
258 jahn 1.13
259     #ifdef CHECK_CONS
260     coj find unused units for darwin_cons output
261     CALL MDSFINDUNIT( DAR_cons_unitP, mythid )
262     open(DAR_cons_unitP,file='darwin_cons_P.txt',status='unknown')
263     CALL MDSFINDUNIT( DAR_cons_unitN, mythid )
264     open(DAR_cons_unitN,file='darwin_cons_N.txt',status='unknown')
265     CALL MDSFINDUNIT( DAR_cons_unitF, mythid )
266     open(DAR_cons_unitF,file='darwin_cons_Fe.txt',status='unknown')
267     CALL MDSFINDUNIT( DAR_cons_unitS, mythid )
268     open(DAR_cons_unitS,file='darwin_cons_Si.txt',status='unknown')
269     #ifdef ALLOW_CARBON
270     CALL MDSFINDUNIT( DAR_cons_unitC, mythid )
271     open(DAR_cons_unitC,file='darwin_cons_C.txt',status='unknown')
272     CALL MDSFINDUNIT( DAR_cons_unitA, mythid )
273     open(DAR_cons_unitA,file='darwin_cons_A.txt',status='unknown')
274     CALL MDSFINDUNIT( DAR_cons_unitO, mythid )
275     open(DAR_cons_unitO,file='darwin_cons_O.txt',status='unknown')
276     #endif
277     #endif
278    
279 jahn 1.3 c myProcId and myThid
280     ENDIF
281 jahn 1.1
282 jahn 1.13 WRITE(msgBuf,'(A)')
283     & '// ======================================================='
284     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
285     & SQUEEZE_RIGHT, myThid )
286     WRITE(msgBuf,'(A)') '// Darwin init variables >>> END <<<'
287     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
288     & SQUEEZE_RIGHT, myThid )
289     WRITE(msgBuf,'(A)')
290     & '// ======================================================='
291     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
292     & SQUEEZE_RIGHT, myThid )
293    
294     _END_MASTER( myThid )
295     _BARRIER
296    
297     C ----------------------------------------------------------------------
298     C From here on we do actual 3d stuff appropriate for init_varia
299    
300     c reduce amount of diaz
301     #ifdef ALLOW_DIAZ
302     IF (nIter0.EQ.PTRACERS_Iter0) THEN
303     do np = 1, npmax
304     if (diazotroph(np) .eq. 1. _d 0) then
305     DO bj = myByLo(myThid), myByHi(myThid)
306     DO bi = myBxLo(myThid), myBxHi(myThid)
307     DO j=1-Oly,sNy+Oly
308     DO i=1-Olx,sNx+Olx
309     DO k=1,nR
310     Ptracer(i,j,k,bi,bj,iPhy+np-1) =
311     & Ptracer(i,j,k,bi,bj,iPhy+np-1)/10. _d 0
312     ENDDO
313     ENDDO
314     ENDDO
315     ENDDO
316     ENDDO
317     endif
318     enddo
319     ENDIF
320     #endif
321    
322     #ifdef GEIDER
323     C this initializes fields...
324     call MONOD_CHECK_CHL(myThid)
325     #endif
326    
327 jahn 1.1 CALL LEF_ZERO( fice,myThid )
328     CALL LEF_ZERO( inputFe,myThid )
329     CALL LEF_ZERO( sur_par,myThid )
330     #ifdef NUT_SUPPLY
331     DO bj = myByLo(myThid), myByHi(myThid)
332     DO bi = myBxLo(myThid), myBxHi(myThid)
333     DO j=1-Oly,sNy+Oly
334     DO i=1-Olx,sNx+Olx
335     DO k=1,nR
336     nut_wvel(i,j,k,bi,bj) = 0. _d 0
337     ENDDO
338     ENDDO
339     ENDDO
340     ENDDO
341     ENDDO
342     #endif
343    
344     #ifdef ALLOW_PAR_DAY
345     DO iPAR=1,2
346     DO bj=myByLo(myThid), myByHi(myThid)
347     DO bi=myBxLo(myThid), myBxHi(myThid)
348     DO k=1,nR
349     DO j=1-Oly,sNy+Oly
350     DO i=1-Olx,sNx+Olx
351     PARday(i,j,k,bi,bj,iPAR) = 0. _d 0
352     ENDDO
353     ENDDO
354     ENDDO
355     ENDDO
356     ENDDO
357     ENDDO
358 jahn 1.10 #endif
359    
360 jahn 1.1 IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
361     & .AND. pickupSuff .EQ. ' ') ) THEN
362     COJ should probably initialize from a file when nIter0 .EQ. 0
363     CALL DARWIN_READ_PICKUP( nIter0, myThid )
364     ENDIF
365     c
366     #ifdef ALLOW_TIMEAVE
367     c set arrays to zero if first timestep
368     DO bj = myByLo(myThid), myByHi(myThid)
369     DO bi = myBxLo(myThid), myBxHi(myThid)
370     CALL TIMEAVE_RESET(PARave, Nr, bi, bj, myThid)
371     CALL TIMEAVE_RESET(PPave, Nr, bi, bj, myThid)
372     CALL TIMEAVE_RESET(Chlave, Nr, bi, bj, myThid)
373     CALL TIMEAVE_RESET(Nfixave, Nr, bi, bj, myThid)
374     CALL TIMEAVE_RESET(Denitave, Nr, bi, bj, myThid)
375 stephd 1.7 #ifdef DAR_DIAG_PARW
376     do i=1,tlam
377     CALL TIMEAVE_RESET(PARwave(1-OLx,1-OLy,1,1,1,i),
378     & Nr,bi,bj,myThid)
379     enddo
380     do np=1,npmax
381     CALL TIMEAVE_RESET(chl2cave(1-OLx,1-OLy,1,1,1,np),
382     & Nr,bi,bj,myThid)
383     enddo
384     #endif
385 jahn 1.1 #ifdef DAR_DIAG_ACDOM
386     CALL TIMEAVE_RESET(aCDOMave, Nr, bi, bj, myThid)
387     #endif
388     #ifdef DAR_DIAG_IRR
389     do i=1,tlam
390     CALL TIMEAVE_RESET(Edave(1-OLx,1-OLy,1,1,1,i),
391     & Nr,bi,bj,myThid)
392     CALL TIMEAVE_RESET(Esave(1-OLx,1-OLy,1,1,1,i),
393     & Nr,bi,bj,myThid)
394     CALL TIMEAVE_RESET(Euave(1-OLx,1-OLy,1,1,1,i),
395     & Nr,bi,bj,myThid)
396 jahn 1.6 CALL TIMEAVE_RESET(Estave(1-OLx,1-OLy,1,1,1,i),
397     & Nr,bi,bj,myThid)
398 jahn 1.1 CALL TIMEAVE_RESET(Eutave(1-OLx,1-OLy,1,1,1,i),
399     & Nr,bi,bj,myThid)
400     enddo
401     #endif
402 jahn 1.5 #ifdef DAR_DIAG_IRR_AMPS
403     do i=1,tlam
404 jahn 1.6 CALL TIMEAVE_RESET(amp1ave(1-OLx,1-OLy,1,1,1,i),
405 jahn 1.5 & Nr,bi,bj,myThid)
406 jahn 1.6 CALL TIMEAVE_RESET(amp2ave(1-OLx,1-OLy,1,1,1,i),
407 jahn 1.5 & Nr,bi,bj,myThid)
408     enddo
409     #endif
410 jahn 1.1 #ifdef DAR_DIAG_ABSORP
411     do i=1,tlam
412     CALL TIMEAVE_RESET(aave(1-OLx,1-OLy,1,1,1,i),
413     & Nr,bi,bj,myThid)
414     enddo
415     #endif
416     #ifdef DAR_DIAG_SCATTER
417     do i=1,tlam
418     CALL TIMEAVE_RESET(btave(1-OLx,1-OLy,1,1,1,i),
419     & Nr,bi,bj,myThid)
420     CALL TIMEAVE_RESET(bbave(1-OLx,1-OLy,1,1,1,i),
421     & Nr,bi,bj,myThid)
422     enddo
423     #endif
424     #ifdef DAR_DIAG_PART_SCATTER
425     do i=1,tlam
426     CALL TIMEAVE_RESET(apartave(1-OLx,1-OLy,1,1,1,i),
427     & Nr,bi,bj,myThid)
428     CALL TIMEAVE_RESET(btpartave(1-OLx,1-OLy,1,1,1,i),
429     & Nr,bi,bj,myThid)
430     CALL TIMEAVE_RESET(bbpartave(1-OLx,1-OLy,1,1,1,i),
431     & Nr,bi,bj,myThid)
432     enddo
433     #endif
434 jahn 1.5 #ifdef DAR_RADTRANS
435     CALL TIMEAVE_RESET(rmudave(1-OLx,1-OLy,1,1),
436     & 1,bi,bj,myThid)
437     #endif
438 jahn 1.1 c ANNA_TAVE
439     #ifdef WAVES_DIAG_PCHL
440     do np=1,npmax
441     CALL TIMEAVE_RESET(Pchlave(1-OLx,1-OLy,1,1,1,np),
442     & Nr,bi,bj,myThid)
443     enddo
444     #endif
445     c ANNA end TAVE
446 stephd 1.7 #ifdef DAR_DIAG_EK
447     do np=1,npmax
448     CALL TIMEAVE_RESET(Ekave(1-OLx,1-OLy,1,1,1,np),
449     & Nr,bi,bj,myThid)
450     CALL TIMEAVE_RESET(EkoverEave(1-OLx,1-OLy,1,1,1,np),
451     & Nr,bi,bj,myThid)
452     do i=1,tlam
453     CALL TIMEAVE_RESET(Ek_nlave(1-OLx,1-OLy,1,1,1,np,i),
454     & Nr,bi,bj,myThid)
455     CALL TIMEAVE_RESET(EkoverE_nlave(1-OLx,1-OLy,1,1,1,np,i),
456     & Nr,bi,bj,myThid)
457     enddo
458     enddo
459     #endif
460 jahn 1.1 #ifdef DAR_DIAG_RSTAR
461     do np=1,npmax
462     CALL TIMEAVE_RESET(Rstarave(1-OLx,1-OLy,1,1,1,np),
463     & Nr,bi,bj,myThid)
464     CALL TIMEAVE_RESET(RNstarave(1-OLx,1-OLy,1,1,1,np),
465     & Nr,bi,bj,myThid)
466     enddo
467     #endif
468     #ifdef DAR_DIAG_DIVER
469     CALL TIMEAVE_RESET(Diver1ave, Nr, bi, bj, myThid)
470     CALL TIMEAVE_RESET(Diver2ave, Nr, bi, bj, myThid)
471     CALL TIMEAVE_RESET(Diver3ave, Nr, bi, bj, myThid)
472     CALL TIMEAVE_RESET(Diver4ave, Nr, bi, bj, myThid)
473     #endif
474     #ifdef ALLOW_DIAZ
475     #ifdef DAR_DIAG_NFIXP
476     do np=1,npmax
477     CALL TIMEAVE_RESET(NfixPave(1-OLx,1-OLy,1,1,1,np),
478     & Nr,bi,bj,myThid)
479     enddo
480     #endif
481     #endif
482     c CALL TIMEAVE_RESET(SURave, 1, bi, bj, myThid)
483    
484     do k=1,Nr
485     DAR_TimeAve(bi,bj,k)=0. _d 0
486     enddo
487     ENDDO
488     ENDDO
489     #endif /* ALLOW_TIMEAVE */
490    
491    
492     RETURN
493     END
494     #endif /*MONOD*/
495     #endif /*ALLOW_PTRACERS*/
496     c ==========================================================
497    

  ViewVC Help
Powered by ViewVC 1.1.22