/[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.3 - (show annotations) (download)
Thu Sep 22 16:01:30 2011 UTC (13 years, 10 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 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 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 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 #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 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 #endif
182 call MONOD_CHECK_CHL(myThid)
183 #endif
184 #endif
185 c ANNA endif
186
187 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
188 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 c myProcId and myThid
283 ENDIF
284
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 IF ( myProcId.EQ.0 .AND. myThid.EQ.1 ) THEN
416 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 c myProcId and myThid
434 ENDIF
435 #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