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

Contents 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.6 - (show annotations) (download)
Thu Aug 23 21:48:25 2012 UTC (12 years, 11 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt64_20121012
Changes since 1.5: +5 -3 lines
add Edstop timeave diagnostic, rename c1/2 to amp1/2

1 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_init_vari.F,v 1.5 2012/07/30 15:21:51 jahn Exp $
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 _BEGIN_MASTER( myThid )
92
93 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 _END_MASTER( myThid )
106
107 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 _BEGIN_MASTER( myThid )
128
129 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 #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 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 #endif
188 #endif
189 #endif
190 _END_MASTER( myThid )
191
192 #ifdef DYNAMIC_CHL
193 C this initializes fields...
194 call MONOD_CHECK_CHL(myThid)
195 #endif
196
197 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
198 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 c myProcId and myThid
293 ENDIF
294
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(Estave(1-OLx,1-OLy,1,1,1,i),
354 & Nr,bi,bj,myThid)
355 CALL TIMEAVE_RESET(Eutave(1-OLx,1-OLy,1,1,1,i),
356 & Nr,bi,bj,myThid)
357 enddo
358 #endif
359 #ifdef DAR_DIAG_IRR_AMPS
360 do i=1,tlam
361 CALL TIMEAVE_RESET(amp1ave(1-OLx,1-OLy,1,1,1,i),
362 & Nr,bi,bj,myThid)
363 CALL TIMEAVE_RESET(amp2ave(1-OLx,1-OLy,1,1,1,i),
364 & Nr,bi,bj,myThid)
365 enddo
366 #endif
367 #ifdef DAR_DIAG_ABSORP
368 do i=1,tlam
369 CALL TIMEAVE_RESET(aave(1-OLx,1-OLy,1,1,1,i),
370 & Nr,bi,bj,myThid)
371 enddo
372 #endif
373 #ifdef DAR_DIAG_SCATTER
374 do i=1,tlam
375 CALL TIMEAVE_RESET(btave(1-OLx,1-OLy,1,1,1,i),
376 & Nr,bi,bj,myThid)
377 CALL TIMEAVE_RESET(bbave(1-OLx,1-OLy,1,1,1,i),
378 & Nr,bi,bj,myThid)
379 enddo
380 #endif
381 #ifdef DAR_DIAG_PART_SCATTER
382 do i=1,tlam
383 CALL TIMEAVE_RESET(apartave(1-OLx,1-OLy,1,1,1,i),
384 & Nr,bi,bj,myThid)
385 CALL TIMEAVE_RESET(btpartave(1-OLx,1-OLy,1,1,1,i),
386 & Nr,bi,bj,myThid)
387 CALL TIMEAVE_RESET(bbpartave(1-OLx,1-OLy,1,1,1,i),
388 & Nr,bi,bj,myThid)
389 enddo
390 #endif
391 #ifdef DAR_RADTRANS
392 CALL TIMEAVE_RESET(rmudave(1-OLx,1-OLy,1,1),
393 & 1,bi,bj,myThid)
394 #endif
395 c ANNA_TAVE
396 #ifdef WAVES_DIAG_PCHL
397 do np=1,npmax
398 CALL TIMEAVE_RESET(Pchlave(1-OLx,1-OLy,1,1,1,np),
399 & Nr,bi,bj,myThid)
400 enddo
401 #endif
402 c ANNA end TAVE
403 #ifdef DAR_DIAG_RSTAR
404 do np=1,npmax
405 CALL TIMEAVE_RESET(Rstarave(1-OLx,1-OLy,1,1,1,np),
406 & Nr,bi,bj,myThid)
407 CALL TIMEAVE_RESET(RNstarave(1-OLx,1-OLy,1,1,1,np),
408 & Nr,bi,bj,myThid)
409 enddo
410 #endif
411 #ifdef DAR_DIAG_DIVER
412 CALL TIMEAVE_RESET(Diver1ave, Nr, bi, bj, myThid)
413 CALL TIMEAVE_RESET(Diver2ave, Nr, bi, bj, myThid)
414 CALL TIMEAVE_RESET(Diver3ave, Nr, bi, bj, myThid)
415 CALL TIMEAVE_RESET(Diver4ave, Nr, bi, bj, myThid)
416 #endif
417 #ifdef ALLOW_DIAZ
418 #ifdef DAR_DIAG_NFIXP
419 do np=1,npmax
420 CALL TIMEAVE_RESET(NfixPave(1-OLx,1-OLy,1,1,1,np),
421 & Nr,bi,bj,myThid)
422 enddo
423 #endif
424 #endif
425 c CALL TIMEAVE_RESET(SURave, 1, bi, bj, myThid)
426 WRITE(msgbuf,'(A)')
427 & 'QQ start timeave'
428 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
429 & SQUEEZE_RIGHT , mythid)
430
431 do k=1,Nr
432 DAR_TimeAve(bi,bj,k)=0. _d 0
433 enddo
434 ENDDO
435 ENDDO
436 #endif /* ALLOW_TIMEAVE */
437
438 #ifdef CHECK_CONS
439 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
440 coj find unused units for darwin_cons output
441 CALL MDSFINDUNIT( DAR_cons_unit1, mythid )
442 open(DAR_cons_unit1,file='darwin_cons_P.txt',status='unknown')
443 CALL MDSFINDUNIT( DAR_cons_unit2, mythid )
444 open(DAR_cons_unit2,file='darwin_cons_N.txt',status='unknown')
445 CALL MDSFINDUNIT( DAR_cons_unit3, mythid )
446 open(DAR_cons_unit3,file='darwin_cons_Fe.txt',status='unknown')
447 CALL MDSFINDUNIT( DAR_cons_unit4, mythid )
448 open(DAR_cons_unit4,file='darwin_cons_Si.txt',status='unknown')
449 #ifdef ALLOW_CARBON
450 CALL MDSFINDUNIT( DAR_cons_unit5, mythid )
451 open(DAR_cons_unit5,file='darwin_cons_C.txt',status='unknown')
452 CALL MDSFINDUNIT( DAR_cons_unit6, mythid )
453 open(DAR_cons_unit6,file='darwin_cons_A.txt',status='unknown')
454 CALL MDSFINDUNIT( DAR_cons_unit7, mythid )
455 open(DAR_cons_unit7,file='darwin_cons_O.txt',status='unknown')
456 #endif
457 c myProcId and myThid
458 ENDIF
459 #endif
460
461 c test....................
462 c write(6,*)'finishing darwin_init_vari '
463 c test....................
464 WRITE(msgBuf,'(A)')
465 & '// ======================================================='
466 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
467 & SQUEEZE_RIGHT, myThid )
468 WRITE(msgBuf,'(A)') '// Darwin init variables >>> END <<<'
469 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
470 & SQUEEZE_RIGHT, myThid )
471 WRITE(msgBuf,'(A)')
472 & '// ======================================================='
473 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
474 & SQUEEZE_RIGHT, myThid )
475
476
477 RETURN
478 END
479 #endif /*MONOD*/
480 #endif /*ALLOW_PTRACERS*/
481 c ==========================================================
482

  ViewVC Help
Powered by ViewVC 1.1.22