/[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.12 - (show annotations) (download)
Fri Jun 21 12:49:47 2013 UTC (12 years, 1 month ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64k_20130723, ctrb_darwin2_ckpt64m_20130820, ctrb_darwin2_ckpt64n_20130826, ctrb_darwin2_ckpt64o_20131024, ctrb_darwin2_ckpt64i_20130622, ctrb_darwin2_ckpt64l_20130806, ctrb_darwin2_ckpt64j_20130704, ctrb_darwin2_ckpt64p_20131118, ctrb_darwin2_ckpt64q_20131118, ctrb_darwin2_ckpt64p_20131024
Changes since 1.11: +15 -15 lines
fix CHECK_CONS io units

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

  ViewVC Help
Powered by ViewVC 1.1.22