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

Annotation of /MITgcm_contrib/darwin2/pkg/monod/monod_forcing.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.17 - (hide annotations) (download)
Fri Sep 12 16:19:39 2014 UTC (10 years, 10 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65j_20150225, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65e_20140929, ctrb_darwin2_ckpt65l_20150504, ctrb_darwin2_ckpt65d_20140915, ctrb_darwin2_ckpt65g_20141120, ctrb_darwin2_ckpt65k_20150402, ctrb_darwin2_ckpt65f_20141014, ctrb_darwin2_ckpt65i_20150123, ctrb_darwin2_ckpt65h_20141217
Changes since 1.16: +48 -2 lines
multiply tendencies by maskInC so tracers are not modified on open boundary

1 jahn 1.17 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/monod/monod_forcing.F,v 1.16 2013/12/27 17:29:00 jahn Exp $
2 jahn 1.2 C $Name: $
3 jahn 1.1
4     #include "CPP_OPTIONS.h"
5     #include "PTRACERS_OPTIONS.h"
6     #include "DARWIN_OPTIONS.h"
7    
8     #ifdef ALLOW_PTRACERS
9     #ifdef ALLOW_MONOD
10    
11     c=============================================================
12     c subroutine MONOD_forcing
13     c step forward bio-chemical tracers in time
14     C==============================================================
15     SUBROUTINE MONOD_FORCING(
16     U Ptr,
17     I bi,bj,imin,imax,jmin,jmax,
18     I myTime,myIter,myThid)
19     #include "SIZE.h"
20     #include "EEPARAMS.h"
21     #include "PARAMS.h"
22     #include "GRID.h"
23     #include "DYNVARS.h"
24 jahn 1.8 c for Qsw and/or surfaceForcingT
25     c choice which field to take pCO2 from for pCO2limit
26     c this assumes we use Ttendency from offline
27 jahn 1.1 #include "FFIELDS.h"
28     #ifdef ALLOW_LONGSTEP
29     #include "LONGSTEP.h"
30     #endif
31     #include "PTRACERS_SIZE.h"
32     #include "PTRACERS_PARAMS.h"
33     #include "GCHEM.h"
34     #include "MONOD_SIZE.h"
35     #include "MONOD.h"
36     #include "DARWIN_IO.h"
37     #include "DARWIN_FLUX.h"
38     #include "MONOD_FIELDS.h"
39    
40     c ANNA include wavebands_params.h
41     #ifdef WAVEBANDS
42     #include "SPECTRAL_SIZE.h"
43     #include "SPECTRAL.h"
44     #include "WAVEBANDS_PARAMS.h"
45     #endif
46    
47 stephd 1.7
48 jahn 1.1 C === Global variables ===
49     c tracers
50     _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)
51     INTEGER bi,bj,imin,imax,jmin,jmax
52     INTEGER myIter
53     _RL myTime
54     INTEGER myThid
55    
56     C !FUNCTIONS:
57     C == Functions ==
58     #ifdef ALLOW_PAR_DAY
59     LOGICAL DIFF_PHASE_MULTIPLE
60     EXTERNAL DIFF_PHASE_MULTIPLE
61     #endif
62    
63     C============== Local variables ============================================
64     c plankton arrays
65     _RL ZooP(nzmax)
66     _RL ZooN(nzmax)
67     _RL ZooFe(nzmax)
68     _RL ZooSi(nzmax)
69     _RL Phy(npmax)
70     _RL Phy_k(npmax,Nr)
71     _RL Phyup(npmax)
72     _RL part_k(Nr)
73 stephd 1.6 #ifdef ALLOW_CDOM
74     _RL cdom_k(Nr)
75     #endif
76 jahn 1.1 c iron partitioning
77     _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78     c some working variables
79     _RL sumpy
80     _RL sumpyup
81     c light variables
82     _RL PAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83     _RL sfac(1-OLy:sNy+OLy)
84     _RL atten,lite
85     _RL newtime ! for sub-timestepping
86     _RL runtim ! time from tracer initialization
87    
88    
89     c ANNA define variables for wavebands
90     #ifdef WAVEBANDS
91     integer ilam
92     _RL PARw_k(tlam,Nr)
93     _RL PARwup(tlam)
94     _RL acdom_k(Nr,tlam)
95 stephd 1.13 _RL Ek_nll(npmax,tlam)
96     _RL EkoverE_nll(npmax,tlam)
97 jahn 1.1 #ifdef DAR_RADTRANS
98     integer iday,iyr,imon,isec,lp,wd,mydate(4)
99     _RL Edwsf(tlam),Eswsf(tlam)
100 jahn 1.10 _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr)
101     _RL Estop(tlam,Nr),Eutop(tlam,Nr)
102 jahn 1.1 _RL tirrq(nr)
103     _RL tirrwq(tlam,nr)
104 jahn 1.10 _RL amp1(tlam,nr), amp2(tlam,nr)
105 jahn 1.1 _RL solz
106     _RL rmud
107     _RL actot,bctot,bbctot
108     _RL apart_k(Nr,tlam),bpart_k(Nr,tlam),bbpart_k(Nr,tlam)
109     _RL bt_k(Nr,tlam), bb_k(Nr,tlam)
110 jahn 1.11 _RL discEs, discEu
111 jahn 1.12 INTEGER idiscEs,jdiscEs,kdiscEs,ldiscEs
112     INTEGER idiscEu,jdiscEu,kdiscEu,ldiscEu
113 jahn 1.1 #else
114     _RL PARwdn(tlam)
115     #endif
116     C always need for diagnostics
117     _RL a_k(Nr,tlam)
118     #endif /* WAVEBANDS */
119    
120    
121     #ifdef DAR_DIAG_DIVER
122     _RL Diver1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
123     _RL Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
124     _RL Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
125     _RL Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
126 jahn 1.5 _RL Shannon(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127     _RL Simpson(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128 jahn 1.1
129     _RL tmpphy(npmax)
130     _RL totphy, biotot, maxphy, phymax
131     #endif
132    
133     #ifdef GEIDER
134     _RL phychl(npmax)
135     _RL phychl_k(npmax,Nr)
136 stephd 1.13 _RL Ekl(npmax)
137     _RL EkoverEl(npmax)
138     _RL chl2cl(npmax)
139 jahn 1.1 #ifdef DYNAMIC_CHL
140     _RL dphychl(npmax)
141     _RL chlup(npmax)
142 stephd 1.13 _RL accliml(npmax)
143 jahn 1.1 #endif
144     #endif
145 stephd 1.6 #ifdef ALLOW_CDOM
146     _RL cdoml
147     _RL dcdoml
148     #endif
149 jahn 1.1
150     #ifdef ALLOW_DIAGNOSTICS
151     COJ for diagnostics
152     _RL PParr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153     _RL Nfixarr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
154     c ANNA_TAVE
155     #ifdef WAVES_DIAG_PCHL
156     _RL Pchlarr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,npmax)
157     #endif
158     c ANNA end TAVE
159     #ifdef DAR_DIAG_RSTAR
160     _RL Rstararr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,npmax)
161     #endif
162     #ifdef ALLOW_DIAZ
163     #ifdef DAR_DIAG_NFIXP
164     _RL NfixParr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,npmax)
165     #endif
166     #endif
167     #endif
168    
169    
170     _RL totphyC
171     #ifdef ALLOW_PAR_DAY
172     LOGICAL itistime
173     INTEGER PARiprev, PARiaccum, iperiod, nav
174     _RL phase
175     _RL dtsubtime
176     #endif
177     #ifdef DAR_DIAG_CHL
178     _RL ChlGeiderlocal, ChlDoneylocal, ChlCloernlocal
179     #ifdef ALLOW_DIAGNOSTICS
180     _RL GeiderChlarr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
181     _RL GeiderChl2Carr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
182     _RL DoneyChlarr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
183     _RL DoneyChl2Carr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
184     _RL CloernChlarr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
185     _RL CloernChl2Carr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
186     #endif
187     #endif
188     c
189     _RL freefu
190     _RL inputFel
191    
192     c some local variables
193     _RL PO4l
194     _RL NO3l
195     _RL FeTl
196     _RL Sil
197     _RL DOPl
198     _RL DONl
199     _RL DOFel
200     _RL POPl
201     _RL PONl
202     _RL POFel
203     _RL PSil
204     _RL POPupl
205     _RL PONupl
206     _RL POFeupl
207     _RL PSiupl
208     _RL Tlocal
209     _RL Slocal
210 stephd 1.7 _RL pCO2local
211 jahn 1.1 _RL Qswlocal
212     _RL NH4l
213     _RL NO2l
214     _RL PARl
215     _RL dzlocal
216     _RL dz_k(Nr)
217     _RL dtplankton
218     _RL bottom
219     _RL PP
220     _RL Nfix
221     _RL denit
222     _RL Chl
223     _RL Rstarl(npmax)
224     _RL RNstarl(npmax)
225     #ifdef DAR_DIAG_GROW
226     _RL Growl(npmax)
227     _RL Growsql(npmax)
228     #endif
229     #ifdef ALLOW_DIAZ
230     #ifdef DAR_DIAG_NFIXP
231     _RL NfixPl(npmax)
232     #endif
233     #endif
234    
235     c local tendencies
236     _RL dphy(npmax)
237     _RL dzoop(nzmax)
238     _RL dzoon(nzmax)
239     _RL dzoofe(nzmax)
240     _RL dzoosi(nzmax)
241     _RL dPO4l
242     _RL dNO3l
243     _RL dFeTl
244     _RL dSil
245     _RL dDOPl
246     _RL dDONl
247     _RL dDOFel
248     _RL dPOPl
249     _RL dPONl
250     _RL dPOFel
251     _RL dPSil
252     _RL dNH4l
253     _RL dNO2l
254    
255     #ifdef ALLOW_CARBON
256     _RL dicl
257     _RL docl
258     _RL pocl
259     _RL picl
260     _RL alkl
261     _RL o2l
262     _RL ZooCl(nzmax)
263     _RL pocupl
264     _RL picupl
265     c tendencies
266     _RL ddicl
267     _RL ddocl
268     _RL dpocl
269     _RL dpicl
270     _RL dalkl
271     _RL do2l
272     _RL dZooCl(nzmax)
273     c air-sea fluxes
274 jahn 1.2 _RL flxCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
275     _RL flxALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
276     _RL flxO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
277 jahn 1.1 #endif
278    
279     _RL tot_Nfix
280    
281     _RL tmp
282    
283     _RL phytmp, chltmp
284    
285     INTEGER i,j,k,it, ktmp
286     INTEGER np, nz, np2, npsave
287     INTEGER debug
288     CHARACTER*8 diagname
289    
290     c
291     c
292     c
293     DO j=1-OLy,sNy+OLy
294     DO i=1-OLx,sNx+OLx
295     do k=1,Nr
296     freefe(i,j,k)=0. _d 0
297     PAR(i,j,k) = 0. _d 0
298     #ifdef DAR_DIAG_DIVER
299     Diver1(i,j,k)=0. _d 0
300     Diver2(i,j,k)=0. _d 0
301     Diver3(i,j,k)=0. _d 0
302     Diver4(i,j,k)=0. _d 0
303 jahn 1.5 Shannon(i,j,k)=0. _d 0
304     Simpson(i,j,k)=1. _d 0
305 jahn 1.1 #endif
306    
307     #ifdef ALLOW_DIAGNOSTICS
308     COJ for diagnostics
309     PParr(i,j,k) = 0. _d 0
310     Nfixarr(i,j,k) = 0. _d 0
311 jahn 1.4 #ifdef DAR_DIAG_CHL
312     GeiderChlarr(i,j,k) = 0. _d 0
313     GeiderChl2Carr(i,j,k) = 0. _d 0
314     DoneyChlarr(i,j,k) = 0. _d 0
315     DoneyChl2Carr(i,j,k) = 0. _d 0
316     CloernChlarr(i,j,k) = 0. _d 0
317     CloernChl2Carr(i,j,k) = 0. _d 0
318     #endif
319 jahn 1.1 c ANNA_TAVE
320     #ifdef WAVES_DIAG_PCHL
321     DO np=1,npmax
322     Pchlarr(i,j,k,np) = 0. _d 0
323     ENDDO
324     #endif
325     c ANNA end TAVE
326     #ifdef DAR_DIAG_RSTAR
327     DO np=1,npmax
328     Rstararr(i,j,k,np) = 0. _d 0
329     ENDDO
330     #endif
331     COJ
332     #ifdef ALLOW_DIAZ
333     #ifdef DAR_DIAG_NFIXP
334     DO np=1,npmax
335     NfixParr(i,j,k,np) = 0. _d 0
336     ENDDO
337     #endif
338     #endif
339     #endif
340     enddo
341     ENDDO
342     ENDDO
343 jahn 1.11
344     #ifdef DAR_RADTRANS
345 jahn 1.12 idiscEs = 0
346     jdiscEs = 0
347     kdiscEs = 0
348     ldiscEs = 0
349     idiscEu = 0
350     jdiscEu = 0
351     kdiscEu = 0
352     ldiscEu = 0
353 jahn 1.11 discEs = 0.
354     discEu = 0.
355     #endif
356 jahn 1.1 c
357     c bio-chemical time loop
358     c--------------------------------------------------
359     DO it=1,nsubtime
360     c -------------------------------------------------
361     tot_Nfix=0. _d 0
362     COJ cannot use dfloat because of adjoint
363     COJ division will be double precision anyway because of dTtracerLev
364     newtime=myTime-dTtracerLev(1)+
365     & float(it)*dTtracerLev(1)/float(nsubtime)
366     c print*,'it ',it,newtime,nsubtime,myTime
367     runtim=myTime-float(PTRACERS_Iter0)*dTtracerLev(1)
368    
369     c determine iron partitioning - solve for free iron
370     c ---------------------------
371     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
372     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
373     & myIter, mythid)
374     c --------------------------
375     #ifdef ALLOW_CARBON
376     c air-sea flux and dilution of CO2
377     call dic_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
378     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
379     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
380     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
381     & flxCO2,
382     & bi,bj,imin,imax,jmin,jmax,
383     & myIter,myTime,myThid)
384     c air-sea flux of O2
385     call dic_o2_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iO2),
386     & flxO2,
387     & bi,bj,imin,imax,jmin,jmax,
388     & myIter,myTime,myThid)
389     c dilusion of alkalinity
390     call dic_alk_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
391     & flxALK,
392     & bi,bj,imin,imax,jmin,jmax,
393     & myIter,myTime,myThid)
394     #endif
395    
396    
397     c find light in each grid cell
398     c ---------------------------
399     c determine incident light
400     #ifndef READ_PAR
401     #ifndef USE_QSW
402     DO j=1-OLy,sNy+OLy
403     sfac(j)=0. _d 0
404     ENDDO
405     call darwin_insol(newTime,sfac,bj)
406     #endif /* not USE_QSW */
407     #endif /* not READ_PAR */
408    
409     #ifdef ALLOW_PAR_DAY
410     C find out which slot of PARday has previous day's average
411     dtsubtime = dTtracerLev(1)/float(nsubtime)
412     C running index of averaging period
413     C myTime has already been incremented in this iteration,
414     C go back half a substep to avoid roundoff problems
415     iperiod = FLOOR((newtime-0.5 _d 0*dtsubtime)
416     & /darwin_PARavPeriod)
417     C 0 -> 1, 1->2, 2->0, ...
418     PARiprev = MOD(iperiod, 2) + 1
419    
420     #ifdef ALLOW_DIAGNOSTICS
421     C always fill; this will be the same during PARavPeriod, but this
422     C way it won't blow up for weird diagnostics periods.
423     C we fill before updating, so the diag is the one used in this time
424     C step
425     CALL DIAGNOSTICS_FILL(
426     & PARday(1-Olx,1-Oly,1,bi,bj,PARiprev),'PARday ',
427     & 0,Nr,2,bi,bj,myThid )
428     #endif
429     #endif /* ALLOW_PAR_DAY */
430    
431     #ifdef DAR_RADTRANS
432     #ifndef DAR_RADTRANS_USE_MODEL_CALENDAR
433     #ifdef ALLOW_CAL
434     C get current date and time of day: iyr/imon/iday+isec
435     CALL CAL_GETDATE( myIter, newtime, mydate, mythid )
436     CALL CAL_CONVDATE( mydate,iyr,imon,iday,isec,lp,wd,mythid )
437     #else
438     STOP 'need cal package or DAR_RADTRANS_USE_MODEL_CALENDAR'
439     #endif
440     #endif
441     #endif
442    
443     C.................................................................
444     C.................................................................
445    
446    
447     C ========================== i,j loops =================================
448     DO j=1,sNy
449     DO i=1,sNx
450    
451     c ------------ these are convenient ------------------------------------
452     DO k=1,Nr
453     part_k(k) = max(Ptr(i,j,k,bi,bj,iPOP),0. _d 0)
454 stephd 1.6 #ifdef ALLOW_CDOM
455     cdom_k(k) = max(Ptr(i,j,k,bi,bj,iCDOM),0. _d 0)
456     #endif
457 jahn 1.1 DO np = 1,npmax
458     Phy_k(np,k) = max(Ptr(i,j,k,bi,bj,iPhy+np-1),0. _d 0)
459     #ifdef GEIDER
460     #ifdef DYNAMIC_CHL
461     phychl_k(np,k) = max(Ptr(i,j,k,bi,bj,iChl+np-1),0. _d 0)
462     #else
463     phychl_k(np,k) = max(Chl_phy(i,j,k,bi,bj,np), 0. _d 0)
464     #endif
465     #endif
466     ENDDO
467     ENDDO
468    
469     c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
470     #ifdef WAVEBANDS
471     #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
472 stephd 1.6 #ifdef ALLOW_CDOM
473     call MONOD_ACDOM(cdom_k,
474     O acdom_k,
475     I myThid)
476     #else
477 jahn 1.1 call MONOD_ACDOM(phychl_k,aphy_chl,aw,
478     O acdom_k,
479     I myThid)
480 stephd 1.6 #endif
481 jahn 1.1 #else
482     DO k=1,Nr
483     DO ilam = 1,tlam
484     acdom_k(k,ilam) = acdom(ilam)
485     ENDDO
486     ENDDO
487     #endif /* DAR_CALC_ACDOM or DAR_RADTRANS */
488     #endif /* WAVEBANDS */
489    
490     c ------------ GET INCIDENT NON-SPECTRAL LIGHT -------------------------
491     #if !(defined(WAVEBANDS) && defined(OASIM))
492     #ifdef READ_PAR
493    
494     lite = sur_par(i,j,bi,bj)
495    
496     #else /* not READ_PAR */
497     #ifdef USE_QSW
498    
499     #ifdef ALLOW_LONGSTEP
500     Qswlocal=LS_Qsw(i,j,bi,bj)
501     #else
502     Qswlocal=Qsw(i,j,bi,bj)
503     #endif
504     lite = -parfrac*Qswlocal*parconv*maskC(i,j,1,bi,bj)
505    
506     #else /* not USE_QSW */
507    
508 jahn 1.14 C convert W/m2 to uEin/s/m2
509     lite = sfac(j)*parconv*maskC(i,j,1,bi,bj)
510 jahn 1.1
511     #endif /* not USE_QSW */
512     #endif /* not READ_PAR */
513    
514     c take ice coverage into account
515     c unless already done in seaice package
516     #if !(defined (ALLOW_SEAICE) && defined (USE_QSW))
517     lite = lite*(1. _d 0-fice(i,j,bi,bj))
518     #endif
519     #endif /* not(WAVEBANDS and OASIM) */
520    
521     c ------------ LIGHT ATTENUATION: --------------------------------------
522     #ifndef WAVEBANDS
523     c ------------ SINGLE-BAND ATTENUATION ---------------------------------
524     atten=0. _d 0
525     do k=1,Nr
526     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
527     sumpyup = sumpy
528     sumpy = 0. _d 0
529     do np=1,npmax
530     #ifdef GEIDER
531     sumpy = sumpy + phychl_k(np,k)
532     #else
533     sumpy = sumpy + Phy_k(np,k)
534     #endif
535     enddo
536     atten= atten + (k0 + kc*sumpy)*5. _d -1*drF(k)
537     if (k.gt.1)then
538     atten = atten + (k0+kc*sumpyup)*5. _d -1*drF(k-1)
539     endif
540     PAR(i,j,k) = lite*exp(-atten)
541     endif
542     enddo
543    
544     #else /* WAVEBANDS */
545     #ifndef DAR_RADTRANS
546     c ------------ WAVEBANDS W/O RADTRANS ----------------------------------
547     do ilam = 1,tlam
548     #ifdef OASIM
549     c add direct and diffuse, convert to uEin/m2/s/nm
550     PARwup(ilam) = WtouEins(ilam)*(oasim_ed(i,j,ilam,bi,bj)+
551     & oasim_es(i,j,ilam,bi,bj))
552     c and take ice fraction into account
553     c PARwup(ilam) = PARwup(ilam)*(1 _d 0 - fice(i,j,bi,bj))
554     #else
555     c sf is per nm; convert to per waveband
556     PARwup(ilam) = wb_width(ilam)*sf(ilam)*lite
557     #endif
558     enddo
559    
560     do k=1,Nr
561     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
562     do ilam = 1,tlam
563     sumpy = 0.
564     do np = 1,npmax
565     c get total attenuation (absorption) by phyto at each wavelength
566     sumpy = sumpy + (phychl_k(np,k)*aphy_chl(np,ilam))
567     enddo
568     c for diagnostic
569     a_k(k,ilam) = aw(ilam) + sumpy + acdom_k(k,ilam)
570     atten = a_k(k,ilam)*drF(k)
571     PARwdn(ilam) = PARwup(ilam)*exp(-atten)
572     enddo
573    
574     c find for the midpoint of the gridcell (gridcell mean)
575     do ilam = 1,tlam
576     C PARw_k(ilam,k)=exp((log(PARwup(ilam))+log(PARwdn(ilam)))*0.5)
577     PARw_k(ilam,k)=sqrt(PARwup(ilam)*PARwdn(ilam))
578     enddo
579    
580     c cycle
581     do ilam=1,tlam
582     PARwup(ilam) = PARwdn(ilam)
583     enddo
584     else
585     do ilam=1,tlam
586     PARw_k(ilam,k) = 0. _d 0
587     enddo
588     endif
589    
590     c sum wavebands for total PAR at the mid point of the gridcell (PARl)
591     PAR(i,j,k) = 0.
592     do ilam = 1,tlam
593     PAR(i,j,k) = PAR(i,j,k) + PARw_k(ilam,k)
594     enddo
595     enddo
596    
597     #else /* DAR_RADTRANS */
598     c ------------ FULL RADIATIVE TRANSFER CODE ----------------------------
599     do ilam = 1,tlam
600     Edwsf(ilam) = oasim_ed(i,j,ilam,bi,bj)
601     Eswsf(ilam) = oasim_es(i,j,ilam,bi,bj)
602     enddo
603    
604     #ifdef DAR_RADTRANS_USE_MODEL_CALENDAR
605     C simplified solar zenith angle for 360-day year and daily averaged light
606     C cos(solz) is average over daylight period
607     call darwin_solz360(newtime, YC(i,j,bi,bj),
608     O solz)
609    
610     #else /* not DAR_RADTRANS_USE_MODEL_CALENDAR */
611     C use calendar date for full solar zenith angle computation
612     C oj: average light effective at noon?
613     solz = 0.0 _d 0
614     isec = 12*3600
615     call radtrans_sfcsolz(rad,iyr,imon,iday,isec,
616     I XC(i,j,bi,bj),YC(i,j,bi,bj),
617     O solz)
618     #endif /* not DAR_RADTRANS_USE_MODEL_CALENDAR */
619    
620     c have Ed,Es below surface - no need for this adjustment on Ed Es for surface affects
621     c do ilam=1,tlam
622     c rod(ilam) = 0.0 _d 0
623     c ros(ilam) = 0.0 _d 0
624     c enddo
625    
626     c compute 1/cos(zenith) for direct light below surface
627     call radtrans_sfcrmud(rad,solz,
628     O rmud)
629    
630     C compute absorption/scattering coefficients for radtrans
631     DO k=1,Nr
632     dz_k(k) = drF(k)*HFacC(i,j,k,bi,bj)
633     DO ilam = 1,tlam
634     c absorption by phyto
635     actot = 0.0
636     bctot = 0.0
637     bbctot = 0.0
638     DO np = 1,npmax
639     actot = actot + phychl_k(np,k)*aphy_chl(np,ilam)
640     bctot = bctot + phychl_k(np,k)*bphy_chl(np,ilam)
641     bbctot = bbctot + phychl_k(np,k)*bbphy_chl(np,ilam)
642     ENDDO
643     c particulate
644     apart_k(k,ilam) = part_k(k)*apart_P(ilam)
645     bpart_k(k,ilam) = part_k(k)*bpart_P(ilam)
646     bbpart_k(k,ilam) = part_k(k)*bbpart_P(ilam)
647     c add water and CDOM
648     a_k(k,ilam) = aw(ilam)+acdom_k(k,ilam)+actot+apart_k(k,ilam)
649     bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
650     bb_k(k,ilam) = darwin_bbw*bw(ilam)+bbctot+bbpart_k(k,ilam)
651     bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
652 jahn 1.10 c initialize output variables
653     Edz(ilam,k) = 0.0
654     Esz(ilam,k) = 0.0
655     Euz(ilam,k) = 0.0
656     Estop(ilam,k) = 0.0
657     Eutop(ilam,k) = 0.0
658     amp1(ilam,k) = 0.0
659     amp2(ilam,k) = 0.0
660 jahn 1.1 ENDDO
661     ENDDO
662    
663 jahn 1.11 IF (darwin_radtrans_niter.GE.0) THEN
664     call MONOD_RADTRANS_ITER(
665 jahn 1.1 I dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
666     I darwin_radtrans_kmax,darwin_radtrans_niter,
667     O Edz,Esz,Euz,Eutop,
668     O tirrq,tirrwq,
669 jahn 1.10 O amp1,amp2,
670 jahn 1.1 I myThid)
671 jahn 1.11 ELSEIF (darwin_radtrans_niter.EQ.-1) THEN
672 jahn 1.1 c dzlocal ?????
673 jahn 1.11 call MONOD_RADTRANS(
674 jahn 1.1 I drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
675     O Edz,Esz,Euz,Eutop,
676     O tirrq,tirrwq,
677     I myThid)
678 jahn 1.11 ELSE
679     call MONOD_RADTRANS_DIRECT(
680     I dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
681     I darwin_radtrans_kmax,
682     O Edz,Esz,Euz,Estop,Eutop,
683     O tirrq,tirrwq,
684     O amp1,amp2,
685     I myThid)
686     #ifdef DAR_CHECK_IRR_CONT
687 jahn 1.12 IF( dz_k(1) .GT. 0.0 )THEN
688 jahn 1.11 DO ilam = 1,tlam
689 jahn 1.12 IF(Eswsf(ilam).GE.darwin_radmodThresh .OR.
690     & Edwsf(ilam).GE.darwin_radmodThresh ) THEN
691     IF(ABS(Estop(ilam,1)-Eswsf(ilam)) .GT. discEs )THEN
692     discEs = ABS(Estop(ilam,1)-Eswsf(ilam))
693     idiscEs = i
694     jdiscEs = j
695     kdiscEs = 1
696     ldiscEs = ilam
697     ENDIF
698     DO k=1,darwin_radtrans_kmax-1
699     IF(ABS(Estop(ilam,k+1)-Esz(ilam,k)) .GT. discEs)THEN
700     discEs = ABS(Estop(ilam,k+1)-Esz(ilam,k))
701     idiscEs = i
702     jdiscEs = j
703     kdiscEs = k+1
704     ldiscEs = ilam
705     ENDIF
706     IF(ABS(Eutop(ilam,k+1)-Euz(ilam,k)) .GT. discEu)THEN
707     discEu = ABS(Eutop(ilam,k+1)-Euz(ilam,k))
708     idiscEu = i
709     jdiscEu = j
710     kdiscEu = k+1
711     ldiscEu = ilam
712     ENDIF
713 jahn 1.11 ENDDO
714 jahn 1.12 ENDIF
715 jahn 1.11 ENDDO
716 jahn 1.12 ENDIF
717     #endif
718 jahn 1.11 ENDIF
719 jahn 1.1 c
720     c uses chl from prev timestep (as wavebands does)
721     c keep like this in case need to consider upwelling irradiance as affecting the grid box above
722     c will pass to plankton: PARw only, but will be for this timestep for RT and prev timestep for WAVBANDS
723     c
724     c now copy
725     DO k=1,Nr
726     PAR(i,j,k) = tirrq(k)
727     DO ilam = 1,tlam
728     PARw_k(ilam,k) = tirrwq(ilam,k)
729     ENDDO
730     ENDDO
731     #endif /* DAR_RADTRANS */
732    
733     c oj: ???
734     c so PARw and PARwup from WAVEBANDS_1D are from previous timestep (attenuation done in plankton)
735     c but PARw and PARwup from WAVEBANDS_3D and RADTRANS are for the current timestep
736    
737     #endif /* WAVEBANDS */
738    
739     C ============================ k loop ==================================
740     c for each layer ...
741     do k= 1, NR
742     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
743    
744     c make sure we only deal with positive definite numbers
745     c brute force...
746     po4l = max(Ptr(i,j,k,bi,bj,iPO4 ),0. _d 0)
747     no3l = max(Ptr(i,j,k,bi,bj,iNO3 ),0. _d 0)
748     fetl = max(Ptr(i,j,k,bi,bj,iFeT ),0. _d 0)
749     sil = max(Ptr(i,j,k,bi,bj,iSi ),0. _d 0)
750     dopl = max(Ptr(i,j,k,bi,bj,iDOP ),0. _d 0)
751     donl = max(Ptr(i,j,k,bi,bj,iDON ),0. _d 0)
752     dofel = max(Ptr(i,j,k,bi,bj,iDOFe ),0. _d 0)
753     DO nz = 1,nzmax
754     ZooP(nz) = max(Ptr(i,j,k,bi,bj,iZooP (nz)),0. _d 0)
755     ZooN(nz) = max(Ptr(i,j,k,bi,bj,iZooN (nz)),0. _d 0)
756     ZooFe(nz) = max(Ptr(i,j,k,bi,bj,iZooFe(nz)),0. _d 0)
757     ZooSi(nz) = max(Ptr(i,j,k,bi,bj,iZooSi(nz)),0. _d 0)
758     ENDDO
759     popl = max(Ptr(i,j,k,bi,bj,iPOP ),0. _d 0)
760     ponl = max(Ptr(i,j,k,bi,bj,iPON ),0. _d 0)
761     pofel = max(Ptr(i,j,k,bi,bj,iPOFe ),0. _d 0)
762     psil = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
763     NH4l = max(Ptr(i,j,k,bi,bj,iNH4 ),0. _d 0)
764     NO2l = max(Ptr(i,j,k,bi,bj,iNO2 ),0. _d 0)
765 stephd 1.6 #ifdef ALLOW_CDOM
766     cdoml = max(Ptr(i,j,k,bi,bj,iCDOM ),0. _d 0)
767     #endif
768 jahn 1.1 #ifdef ALLOW_CARBON
769     dicl = max(Ptr(i,j,k,bi,bj,iDIC ),0. _d 0)
770     docl = max(Ptr(i,j,k,bi,bj,iDOC ),0. _d 0)
771     pocl = max(Ptr(i,j,k,bi,bj,iPOC ),0. _d 0)
772     picl = max(Ptr(i,j,k,bi,bj,iPIC ),0. _d 0)
773     alkl = max(Ptr(i,j,k,bi,bj,iALK ),0. _d 0)
774     o2l = max(Ptr(i,j,k,bi,bj,iO2 ),0. _d 0)
775     DO nz = 1,nzmax
776     ZooCl(nz) = max(Ptr(i,j,k,bi,bj,iZooC (nz)),0. _d 0)
777     ENDDO
778     #endif
779    
780     totphyC = 0. _d 0
781     DO np=1,npmax
782     totphyC = totphyC + R_PC(np)*Ptr(i,j,k,bi,bj,iPhy+np-1)
783     ENDDO
784    
785     DO np = 1,npmax
786     Phy(np) = Phy_k(np,k)
787     #ifdef GEIDER
788     phychl(np) = phychl_k(np,k)
789     #endif
790     ENDDO
791    
792     #ifdef DAR_DIAG_DIVER
793     Diver1(i,j,k)=0. _d 0
794     Diver2(i,j,k)=0. _d 0
795     Diver3(i,j,k)=0. _d 0
796     Diver4(i,j,k)=0. _d 0
797     totphy=0. _d 0
798     do np=1,npmax
799     totphy=totphy + Phy(np)
800     tmpphy(np)=Phy(np)
801     enddo
802     if (totphy.gt.diver_thresh0) then
803     do np=1,npmax
804     c simple threshhold
805     if (Phy(np).gt.diver_thresh1) then
806     Diver1(i,j,k)=Diver1(i,j,k)+1. _d 0
807     endif
808     c proportion of total biomass
809     if (Phy(np)/totphy.gt.diver_thresh2) then
810     Diver2(i,j,k)=Diver2(i,j,k)+1. _d 0
811     endif
812     enddo
813     c majority of biomass by finding rank order
814     biotot=0. _d 0
815     do np2=1,npmax
816     phymax=0. _d 0
817     do np=1,npmax
818     if (tmpphy(np).gt.phymax) then
819     phymax=tmpphy(np)
820     npsave=np
821     endif
822     enddo
823     if (biotot.lt.totphy*diver_thresh3) then
824     Diver3(i,j,k)=Diver3(i,j,k)+1. _d 0
825     endif
826     biotot=biotot+tmpphy(npsave)
827     tmpphy(npsave)=0. _d 0
828     if (np2.eq.1) then
829     maxphy=phymax
830     endif
831     enddo
832     c ratio of maximum species
833     do np=1,npmax
834     if (Phy(np).gt.diver_thresh4*maxphy) then
835     Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
836     endif
837     enddo
838 jahn 1.5 c totphy > thresh0
839     endif
840     c Shannon and Simpson indices
841     Shannon(i,j,k) = 0. _d 0
842     c note: minimal valid value is 1, but we set to zero below threshold
843     Simpson(i,j,k) = 0. _d 0
844     if (totphy.gt.shannon_thresh) then
845     do np=1,npmax
846     if (Phy(np) .gt. 0. _d 0) then
847     tmpphy(np) = Phy(np)/totphy
848     Shannon(i,j,k)=Shannon(i,j,k)+tmpphy(np)*LOG(tmpphy(np))
849     Simpson(i,j,k)=Simpson(i,j,k)+tmpphy(np)*tmpphy(np)
850     endif
851     enddo
852     Shannon(i,j,k) = -Shannon(i,j,k)
853     Simpson(i,j,k) = 1./Simpson(i,j,k)
854 jahn 1.1 endif
855     #endif
856    
857     c..........................................................
858     c find local light
859     c..........................................................
860    
861     PARl = PAR(i,j,k)
862     c..........................................................
863    
864     c for explicit sinking of particulate matter and phytoplankton
865     if (k.eq.1) then
866     popupl =0. _d 0
867     ponupl =0. _d 0
868     pofeupl = 0. _d 0
869     psiupl = 0. _d 0
870     do np=1,npmax
871     Phyup(np)=0. _d 0
872     #ifdef DYNAMIC_CHL
873     chlup(np)=0. _d 0
874     #endif
875     enddo
876     #ifdef ALLOW_CARBON
877     pocupl = 0. _d 0
878     picupl = 0. _d 0
879     #endif
880     endif
881    
882     #ifdef ALLOW_LONGSTEP
883     Tlocal = LS_theta(i,j,k,bi,bj)
884     Slocal = LS_salt(i,j,k,bi,bj)
885     #else
886     Tlocal = theta(i,j,k,bi,bj)
887     Slocal = salt(i,j,k,bi,bj)
888     #endif
889    
890 stephd 1.7 c choice where to get pCO2 from
891     c taking from igsm dic run - fed through Tflux array
892     c pCO2local=surfaceForcingT(i,j,bi,bj)
893     c or from darwin carbon module
894     #ifdef ALLOW_CARBON
895 stephd 1.15 #ifdef pH_3D
896     pCO2local=pCO2(i,j,k,bi,bj)
897     #else
898 stephd 1.7 pCO2local=pCO2(i,j,bi,bj)
899 stephd 1.15 #endif
900 stephd 1.7 #else
901     pCO2local=280. _d -6
902     #endif
903    
904 jahn 1.1 freefu = max(freefe(i,j,k),0. _d 0)
905     if (k.eq.1) then
906     inputFel = inputFe(i,j,bi,bj)
907     else
908     inputFel = 0. _d 0
909     endif
910    
911     dzlocal = drF(k)*HFacC(i,j,k,bi,bj)
912     c set bottom=1.0 if the layer below is not ocean
913     ktmp=min(nR,k+1)
914     if(hFacC(i,j,ktmp,bi,bj).eq.0. _d 0.or.k.eq.Nr) then
915     bottom = 1.0 _d 0
916     else
917     bottom = 0.0 _d 0
918     endif
919    
920     c set tendencies to 0
921     do np=1,npmax
922     dphy(np)=0. _d 0
923     enddo
924     do nz=1,nzmax
925     dzoop(nz)=0. _d 0
926     dzoon(nz)=0. _d 0
927     dzoofe(nz)=0. _d 0
928     dzoosi(nz)=0. _d 0
929     enddo
930     dPO4l=0. _d 0
931     dNO3l=0. _d 0
932     dFeTl=0. _d 0
933     dSil=0. _d 0
934     dDOPl=0. _d 0
935     dDONl=0. _d 0
936     dDOFel=0. _d 0
937     dPOPl=0. _d 0
938     dPONl=0. _d 0
939     dPOFel=0. _d 0
940     dPSil=0. _d 0
941     dNH4l=0. _d 0
942     dNO2l=0. _d 0
943     #ifdef DYNAMIC_CHL
944     do np=1,npmax
945     dphychl(np)=0. _d 0
946     enddo
947     #endif
948 stephd 1.6 #ifdef ALLOW_CDOM
949     dcdoml=0. _d 0
950     #endif
951 jahn 1.1 #ifdef ALLOW_CARBON
952     ddicl=0. _d 0
953     ddocl=0. _d 0
954     dpocl=0. _d 0
955     dpicl=0. _d 0
956     dalkl=0. _d 0
957     do2l=0. _d 0
958     do nz=1,nzmax
959     dzoocl(nz)=0. _d 0
960     enddo
961     #endif
962     c set other arguments to zero
963     PP=0. _d 0
964     Nfix=0. _d 0
965     denit=0. _d 0
966     do np=1,npmax
967     Rstarl(np)=0. _d 0
968     RNstarl(np)=0. _d 0
969     #ifdef DAR_DIAG_GROW
970     Growl(np)=0. _d 0
971     Growsql(np)=0. _d 0
972     #endif
973     #ifdef ALLOW_DIAZ
974     #ifdef DAR_DIAG_NFIXP
975     NfixPl(np)=0. _d 0
976     #endif
977     #endif
978 stephd 1.13 #ifdef DAR_DIAG_PARW
979     chl2cl(np)=0. _d 0
980     #endif
981     #ifdef DAR_DIAG_EK
982     Ekl(np)=0. _d 0
983     EkoverEl(np)=0. _d 0
984     do ilam=1,tlam
985     Ek_nll(np,ilam)=0. _d 0
986     EkoverE_nll(np,ilam)=0. _d 0
987     enddo
988     #endif
989 jahn 1.1 enddo
990    
991    
992     debug=0
993     c if (i.eq.20.and.j.eq.20.and.k.eq.1) debug=8
994     c if (i.eq.10.and.j.eq.10.and.k.eq.1) debug=100
995     c if (i.eq.1.and.j.eq.10.and.k.eq.1) debug=10
996     c if (i.eq.1.and.j.eq.1.and.k.eq.10) debug=14
997    
998     if (debug.eq.7) print*,'PO4, DOP, POP, ZooP',
999     & PO4l, DOPl, POPl, zooP
1000     if (debug.eq.7) print*,'NO3, NO2, NH4, DON, PON, ZooN',
1001     & NO3l,NO2l,NH4l, DONl, PONl, ZooN
1002     if (debug.eq.7) print*,'FeT, DOFe, POFe, Zoofe',
1003     & FeTl, DOFel, POFel, zooFe
1004     if (debug.eq.7) print*,'Si, Psi, zooSi',
1005     & Sil, PSil, zooSi
1006     if (debug.eq.7) print*,'Total Phy', sumpy, PARl, lite
1007     if (debug.eq.7) print*,'Phy', Phy
1008    
1009     if (debug.eq.8) print*,'k, PARl, inputFel, dzlocal',
1010     & PARl, inputFel, dzlocal
1011    
1012     c if (NO3l.eq.0. _d 0.or.NO2l.eq.0. _d 0
1013     c & .or.NH4l.eq.0. _d 0) then
1014     c print*,'QQ N zeros',i,j,k,NO3l,NO2l,NH4l
1015     c endif
1016    
1017    
1018     c ANNA pass extra variables if WAVEBANDS
1019     CALL MONOD_PLANKTON(
1020     U Phy,
1021     I zooP, zooN, zooFe, zooSi,
1022     O PP, Chl, Nfix, denit,
1023     I PO4l, NO3l, FeTl, Sil,
1024     I NO2l, NH4l,
1025     I DOPl, DONl, DOFel,
1026     I POPl, PONl, POFel, PSil,
1027     I phyup, popupl, ponupl,
1028     I pofeupl, psiupl,
1029     I PARl,
1030     I Tlocal, Slocal,
1031 stephd 1.7 I pCO2local,
1032 jahn 1.1 I freefu, inputFel,
1033     I bottom, dzlocal,
1034     O Rstarl, RNstarl,
1035     #ifdef DAR_DIAG_GROW
1036     O Growl, Growsql,
1037     #endif
1038     #ifdef ALLOW_DIAZ
1039     #ifdef DAR_DIAG_NFIXP
1040     O NfixPl,
1041     #endif
1042     #endif
1043     O dphy, dzooP, dzooN, dzooFe,
1044     O dzooSi,
1045     O dPO4l, dNO3l, dFeTl, dSil,
1046     O dNH4l, dNO2l,
1047     O dDOPl, dDONl, dDOFel,
1048     O dPOPl, dPONl, dPOFel, dPSil,
1049     #ifdef ALLOW_CARBON
1050     I dicl, docl, pocl, picl,
1051     I alkl, o2l, zoocl,
1052     I pocupl, picupl,
1053     O ddicl, ddocl, dpocl, dpicl,
1054     O dalkl, do2l, dzoocl,
1055     #endif
1056     #ifdef GEIDER
1057     O phychl,
1058 stephd 1.13 #ifdef DAR_DIAG_EK
1059     I Ekl, EkoverEl,
1060     #endif
1061     #ifdef DAR_DIAG_PARW
1062     I chl2cl,
1063     #endif
1064 jahn 1.1 #ifdef DYNAMIC_CHL
1065     I dphychl,
1066     I chlup,
1067 stephd 1.13 #ifdef DAR_DIAG_EK
1068     O accliml,
1069     #endif
1070 jahn 1.1 #endif
1071 stephd 1.6 #ifdef ALLOW_CDOM
1072     O dcdoml,
1073     I cdoml,
1074     #endif
1075 jahn 1.1 #ifdef WAVEBANDS
1076     I PARw_k(1,k),
1077 stephd 1.13 #ifdef DAR_DIAG_EK
1078     I Ek_nll, EkoverE_nll,
1079     #endif
1080 jahn 1.1 #endif
1081     #endif
1082     #ifdef ALLOW_PAR_DAY
1083     I PARday(i,j,k,bi,bj,PARiprev),
1084     #endif
1085     #ifdef DAR_DIAG_CHL
1086     O ChlGeiderlocal, ChlDoneylocal,
1087     O ChlCloernlocal,
1088     #endif
1089     I debug,
1090     I runtim,
1091     I MyThid)
1092    
1093     c
1094     c if (i.eq.1.and.k.eq.1.and.j.eq.5) then
1095     c print*,i,j,k
1096     c print*,'NO3,No2,NH4', NO3l, NO2l, NH4l
1097     c print*,'dNO3 etc',dNO3l,dNH4l, dNO2l
1098     c print*,'PO4',PO4l,dPO4l
1099     c endif
1100     c
1101     #ifdef IRON_SED_SOURCE
1102     c only above minimum depth (continental shelf)
1103 jahn 1.3 if (rF(k).gt.-depthfesed) then
1104 jahn 1.1 c only if bottom layer
1105     if (bottom.eq.1.0 _d 0) then
1106     #ifdef IRON_SED_SOURCE_VARIABLE
1107     c calculate sink of POP into bottom layer
1108     tmp=(wp_sink*POPupl)/(dzlocal)
1109     c convert to dPOCl
1110     dFetl=dFetl+fesedflux_pcm*(tmp*106. _d 0)
1111     #else
1112     dFetl=dFetl+fesedflux/
1113     & (drF(k)*hFacC(i,j,k,bi,bj))
1114     #endif
1115     endif
1116     endif
1117     #endif
1118    
1119    
1120     popupl = POPl
1121     ponupl = PONl
1122     pofeupl = POFel
1123     psiupl = PSil
1124     do np=1,npmax
1125     Phyup(np) = Phy(np)
1126     #ifdef DYNAMIC_CHL
1127     chlup(np) = phychl(np)
1128     #endif
1129     enddo
1130    
1131    
1132     c
1133     #ifdef ALLOW_CARBON
1134     pocupl = POCl
1135     picupl = PICl
1136     c include surface forcing
1137     if (k.eq.1) then
1138 jahn 1.2 ddicl = ddicl + flxCO2(i,j)
1139     dalkl = dalkl + flxALK(i,j)
1140     do2l = do2l + flxO2(i,j)
1141 jahn 1.1 endif
1142     #endif
1143     c
1144     #ifdef CONS_SUPP
1145     c only works for two layer model
1146     if (k.eq.2) then
1147     dpo4l=0. _d 0
1148     dno3l=0. _d 0
1149     dfetl=0. _d 0
1150     dsil=0. _d 0
1151     endif
1152     #endif
1153     #ifdef RELAX_NUTS
1154     #ifdef DENIT_RELAX
1155     if (rF(k).lt.-depthdenit) then
1156     if (darwin_relaxscale.gt.0. _d 0) then
1157     IF ( darwin_NO3_RelaxFile .NE. ' ' ) THEN
1158     c Fanny's formulation
1159     tmp=(Ptr(i,j,k,bi,bj,iNO3 )-no3_obs(i,j,k,bi,bj))
1160     if (tmp.gt.0. _d 0) then
1161     dno3l=dno3l-(tmp/
1162     & darwin_relaxscale)
1163     denit=tmp/
1164     & darwin_relaxscale
1165     else
1166     denit=0. _d 0
1167     endif
1168     c --- end fanny's formulation
1169     ENDIF
1170     c steph's alternative
1171     c tmp=(Ptr(i,j,k,bi,bj,iNO3 )-
1172     c & 16. _d 0 * Ptr(i,j,k,bi,bj,iPO4 ))
1173     c if (tmp.gt.0. _d 0) then
1174     c dno3l=dno3l-(tmp/
1175     c & darwin_relaxscale)
1176     c denit=tmp/
1177     c & darwin_relaxscale
1178     c else
1179     c denit=0. _d 0
1180     c endif
1181     c ---- end steph's alternative
1182     endif
1183     endif
1184     #else
1185     if (darwin_relaxscale.gt.0. _d 0) then
1186     IF ( darwin_PO4_RelaxFile .NE. ' ' ) THEN
1187     tmp=(Ptr(i,j,k,bi,bj,iPO4 )-po4_obs(i,j,k,bi,bj))
1188     if (tmp.lt.0. _d 0) then
1189     dpo4l=dpo4l-(tmp/
1190     & darwin_relaxscale)
1191     endif
1192     ENDIF
1193     IF ( darwin_NO3_RelaxFile .NE. ' ' ) THEN
1194     tmp=(Ptr(i,j,k,bi,bj,iNO3 )-no3_obs(i,j,k,bi,bj))
1195     if (tmp.lt.0. _d 0) then
1196     dno3l=dno3l-(tmp/
1197     & darwin_relaxscale)
1198     endif
1199     ENDIF
1200     IF ( darwin_Fet_RelaxFile .NE. ' ' ) THEN
1201     tmp=(Ptr(i,j,k,bi,bj,iFeT )-fet_obs(i,j,k,bi,bj))
1202     if (tmp.lt.0. _d 0) then
1203     dfetl=dfetl-(tmp/
1204     & darwin_relaxscale)
1205     endif
1206     ENDIF
1207     IF ( darwin_Si_RelaxFile .NE. ' ' ) THEN
1208     tmp=( Ptr(i,j,k,bi,bj,iSi )-si_obs(i,j,k,bi,bj))
1209     if (tmp.lt.0. _d 0) then
1210     dsil=dsil-(tmp/
1211     & darwin_relaxscale)
1212     endif
1213     ENDIF
1214     endif
1215     #endif
1216     #endif
1217     #ifdef FLUX_NUTS
1218     dpo4l=dpo4l+po4_flx(i,j,k,bi,bj)
1219     dno3l=dno3l+no3_flx(i,j,k,bi,bj)
1220     dfetl=dfetl+fet_flx(i,j,k,bi,bj)
1221     dsil=dsil+si_flx(i,j,k,bi,bj)
1222     #endif
1223 jahn 1.17
1224     #ifdef ALLOW_OBCS
1225     IF (useOBCS) THEN
1226     dpo4l = dpo4l *maskInC(i,j,bi,bj)
1227     dno3l = dno3l *maskInC(i,j,bi,bj)
1228     dfetl = dfetl *maskInC(i,j,bi,bj)
1229     dsil = dsil *maskInC(i,j,bi,bj)
1230     ddopl = ddopl *maskInC(i,j,bi,bj)
1231     ddonl = ddonl *maskInC(i,j,bi,bj)
1232     ddofel = ddofel*maskInC(i,j,bi,bj)
1233     dpopl = dpopl *maskInC(i,j,bi,bj)
1234     dponl = dponl *maskInC(i,j,bi,bj)
1235     dpofel = dpofel*maskInC(i,j,bi,bj)
1236     dpsil = dpsil *maskInC(i,j,bi,bj)
1237     dnh4l = dnh4l *maskInC(i,j,bi,bj)
1238     dno2l = dno2l *maskInC(i,j,bi,bj)
1239     DO nz = 1,nzmax
1240     dzoop (nz) = dzoop (nz)*maskInC(i,j,bi,bj)
1241     dzoon (nz) = dzoon (nz)*maskInC(i,j,bi,bj)
1242     dzoofe(nz) = dzoofe(nz)*maskInC(i,j,bi,bj)
1243     dzoosi(nz) = dzoosi(nz)*maskInC(i,j,bi,bj)
1244     ENDDO
1245     DO np = 1,npmax
1246     dPhy(np) = dPhy(np)*maskInC(i,j,bi,bj)
1247     #ifdef GEIDER
1248     #ifdef DYNAMIC_CHL
1249     dphychl(np) = dphychl(np)*maskInC(i,j,bi,bj)
1250     #endif
1251     #endif
1252     ENDDO
1253     #ifdef ALLOW_CDOM
1254     dcdoml = dcdoml*maskInC(i,j,bi,bj)
1255     #endif
1256     #ifdef ALLOW_CARBON
1257     ddicl = ddicl*maskInC(i,j,bi,bj)
1258     ddocl = ddocl*maskInC(i,j,bi,bj)
1259     dpocl = dpocl*maskInC(i,j,bi,bj)
1260     dpicl = dpicl*maskInC(i,j,bi,bj)
1261     dalkl = dalkl*maskInC(i,j,bi,bj)
1262     do2l = do2l *maskInC(i,j,bi,bj)
1263     DO nz = 1,nzmax
1264     dzoocl(nz) = dzoocl(nz)*maskInC(i,j,bi,bj)
1265     ENDDO
1266     #endif
1267     ENDIF
1268     #endif
1269    
1270 jahn 1.1 c now update main tracer arrays
1271     dtplankton = PTRACERS_dTLev(k)/float(nsubtime)
1272     Ptr(i,j,k,bi,bj,iPO4 ) = Ptr(i,j,k,bi,bj,iPO4) +
1273     & dtplankton*dpo4l
1274     Ptr(i,j,k,bi,bj,iNO3 ) = Ptr(i,j,k,bi,bj,iNO3) +
1275     & dtplankton*dno3l
1276     Ptr(i,j,k,bi,bj,iFeT ) = Ptr(i,j,k,bi,bj,iFeT) +
1277     & dtplankton*dfetl
1278     Ptr(i,j,k,bi,bj,iSi ) = Ptr(i,j,k,bi,bj,iSi ) +
1279     & dtplankton*dsil
1280     Ptr(i,j,k,bi,bj,iDOP ) = Ptr(i,j,k,bi,bj,iDOP) +
1281     & dtplankton*ddopl
1282     Ptr(i,j,k,bi,bj,iDON ) = Ptr(i,j,k,bi,bj,iDON) +
1283     & dtplankton*ddonl
1284     Ptr(i,j,k,bi,bj,iDOFe) = Ptr(i,j,k,bi,bj,iDOFe) +
1285     & dtplankton*ddofel
1286     Ptr(i,j,k,bi,bj,iPOP ) = Ptr(i,j,k,bi,bj,iPOP ) +
1287     & dtplankton*dpopl
1288     Ptr(i,j,k,bi,bj,iPON ) = Ptr(i,j,k,bi,bj,iPON ) +
1289     & dtplankton*dponl
1290     Ptr(i,j,k,bi,bj,iPOFe) = Ptr(i,j,k,bi,bj,iPOFe) +
1291     & dtplankton*dpofel
1292     Ptr(i,j,k,bi,bj,iPOSi) = Ptr(i,j,k,bi,bj,iPOSi) +
1293     & dtplankton*dpsil
1294     Ptr(i,j,k,bi,bj,iNH4 ) = Ptr(i,j,k,bi,bj,iNH4 ) +
1295     & dtplankton*dnh4l
1296     Ptr(i,j,k,bi,bj,iNO2 ) = Ptr(i,j,k,bi,bj,iNO2 ) +
1297     & dtplankton*dno2l
1298     DO nz = 1,nzmax
1299     Ptr(i,j,k,bi,bj,iZooP (nz)) = Ptr(i,j,k,bi,bj,iZooP (nz)) +
1300     & dtplankton*dzoop (nz)
1301     Ptr(i,j,k,bi,bj,iZooN (nz)) = Ptr(i,j,k,bi,bj,iZooN (nz)) +
1302     & dtplankton*dzoon (nz)
1303     Ptr(i,j,k,bi,bj,iZooFe(nz)) = Ptr(i,j,k,bi,bj,iZooFe(nz)) +
1304     & dtplankton*dzoofe(nz)
1305     Ptr(i,j,k,bi,bj,iZooSi(nz)) = Ptr(i,j,k,bi,bj,iZooSi(nz)) +
1306     & dtplankton*dzoosi(nz)
1307     ENDDO
1308     DO np = 1,npmax
1309     Ptr(i,j,k,bi,bj,iPhy+np-1) = Ptr(i,j,k,bi,bj,iPhy+np-1) +
1310     & dtplankton*dPhy(np)
1311     #ifdef GEIDER
1312     #ifdef DYNAMIC_CHL
1313     if (np.eq.1) Chl=0. _d 0
1314     Ptr(i,j,k,bi,bj,iChl+np-1) = Ptr(i,j,k,bi,bj,iChl+np-1) +
1315     & dtplankton*dphychl(np)
1316     c chltmp=Ptr(i,j,k,bi,bj,iChl+np-1)
1317     c phytmp=Ptr(i,j,k,bi,bj,iPhy+np-1)
1318     c Ptr(i,j,k,bi,bj,iChl+np-1)=
1319     c & max(chltmp,phytmp*R_PC(np)*chl2cmin(np))
1320     c if (np.eq.1.and.i.eq.1.and.j.eq.1.and.k.eq.1)
1321     c & print*,chltmp,phytmp,phytmp*R_PC(np)*chl2cmin(np),
1322     c & phytmp*R_PC(np)*chl2cmax(np)
1323     c in darwin_plankton this is stored for previous timestep. Reset here.
1324     Chl=Chl+Ptr(i,j,k,bi,bj,iChl+np-1)
1325     #else
1326     Chl_phy(i,j,k,bi,bj,np)=phychl(np)
1327     #endif
1328     #endif
1329     ENDDO
1330 stephd 1.6 #ifdef ALLOW_CDOM
1331     Ptr(i,j,k,bi,bj,iCDOM ) = Ptr(i,j,k,bi,bj,iCDOM ) +
1332     & dtplankton*dcdoml
1333     #endif
1334 jahn 1.1 #ifdef ALLOW_CARBON
1335     Ptr(i,j,k,bi,bj,iDIC ) = Ptr(i,j,k,bi,bj,iDIC ) +
1336     & dtplankton*ddicl
1337     Ptr(i,j,k,bi,bj,iDOC ) = Ptr(i,j,k,bi,bj,iDOC ) +
1338     & dtplankton*ddocl
1339     Ptr(i,j,k,bi,bj,iPOC ) = Ptr(i,j,k,bi,bj,iPOC ) +
1340     & dtplankton*dpocl
1341     Ptr(i,j,k,bi,bj,iPIC ) = Ptr(i,j,k,bi,bj,iPIC ) +
1342     & dtplankton*dpicl
1343     Ptr(i,j,k,bi,bj,iALK ) = Ptr(i,j,k,bi,bj,iALK ) +
1344     & dtplankton*dalkl
1345     Ptr(i,j,k,bi,bj,iO2 ) = Ptr(i,j,k,bi,bj,iO2 ) +
1346     & dtplankton*do2l
1347     DO nz = 1,nzmax
1348     Ptr(i,j,k,bi,bj,iZooC (nz)) = Ptr(i,j,k,bi,bj,iZooC (nz)) +
1349     & dtplankton*dzoocl (nz)
1350     ENDDO
1351     #endif
1352     c
1353     #ifdef ALLOW_MUTANTS
1354     cQQQQTEST
1355     if (debug.eq.11) then
1356     if (k.lt.8) then
1357     do np=1,60
1358     if(mod(np,4).eq. 1. _d 0)then
1359     np2=np+1
1360     np4=np+3
1361    
1362     Coj: couldn't test this part after change Phynp -> Ptr(...,iPhy+np-1)
1363     Coj: used to be many copies of this:
1364     C if (dPhy(2).gt.dPhy(4).and.dPhy(4).gt.0. _d 0) then
1365     C print*,'QQQ dphy2 > dphy4',i,j,k,Phy2(i,j,k),
1366     C & Phy4(i,j,k), dPhy(2), dPhy(4)
1367     C endif
1368     C if (Phy2(i,j,k).gt.Phy4(i,j,k).and.
1369     C & Phy4(i,j,k).gt.0. _d 0) then
1370     C print*,'QQ phy02 > phy04',i,j,k,Phy2(i,j,k),
1371     C & Phy4(i,j,k), dPhy(2), dPhy(4)
1372     C endif
1373    
1374     if (dPhy(np2).gt.dPhy(np4).and.dPhy(np4).gt.0. _d 0) then
1375     print*,'QQQ dphy',np2,' > dphy',np4,i,j,k,Phy2(i,j,k),
1376     & Ptr(i,j,k,bi,bj,iPhy+np4-1), dPhy(2), dPhy(4)
1377     endif
1378     if (Ptr(i,j,k,bi,bj,iphy+np2-1).gt.Ptr(i,j,k,bi,bj,iPhy+np4-1)
1379     & .and. Ptr(i,j,k,bi,bj,iPhy+np4-1).gt.0. _d 0) then
1380     print*,'QQ phy',np2,' > ',np4,i,j,k,
1381     & Ptr(i,j,k,bi,bj,iPhy+np2-1),
1382     & Ptr(i,j,k,bi,bj,iPhy+np4-1), dPhy(2), dPhy(4)
1383     endif
1384    
1385     endif
1386     enddo ! np
1387     endif ! k
1388     endif
1389     #endif
1390    
1391     #ifdef ALLOW_DIAGNOSTICS
1392     COJ for diagnostics
1393     PParr(i,j,k) = PP
1394     Nfixarr(i,j,k) = Nfix
1395     c ANNA_TAVE
1396     #ifdef WAVES_DIAG_PCHL
1397     DO np = 1,npmax
1398     Pchlarr(i,j,k,np) = phychl(np)
1399     ENDDO
1400     #endif
1401     c ANNA end TAVE
1402     #ifdef DAR_DIAG_RSTAR
1403     DO np = 1,npmax
1404     Rstararr(i,j,k,np) = Rstarl(np)
1405     ENDDO
1406     #endif
1407     #ifdef ALLOW_DIAZ
1408     #ifdef DAR_DIAG_NFIXP
1409     DO np = 1,npmax
1410     NfixParr(i,j,k,np) = NfixPl(np)
1411     ENDDO
1412     #endif
1413     #endif
1414     #ifdef DAR_DIAG_CHL
1415     GeiderChlarr(i,j,k) = ChlGeiderlocal
1416     DoneyChlarr(i,j,k) = ChlDoneylocal
1417     CloernChlarr(i,j,k) = ChlCloernlocal
1418     IF (totphyC .NE. 0. _d 0) THEN
1419     GeiderChl2Carr(i,j,k) = ChlGeiderlocal/totphyC
1420     DoneyChl2Carr(i,j,k) = ChlDoneylocal/totphyC
1421     CloernChl2Carr(i,j,k) = ChlCloernlocal/totphyC
1422     ELSE
1423     GeiderChl2Carr(i,j,k) = 0. _d 0
1424     DoneyChl2Carr(i,j,k) = 0. _d 0
1425     CloernChl2Carr(i,j,k) = 0. _d 0
1426     ENDIF
1427     #endif
1428     COJ
1429     #endif /* ALLOW_DIAGNOSTICS */
1430    
1431     c total fixation (NOTE - STILL NEEDS GLOB SUM)
1432     tot_Nfix=tot_Nfix+
1433     & Nfix*rA(i,j,bi,bj)*rF(k)*hFacC(i,j,k,bi,bj)
1434    
1435     #ifdef ALLOW_TIMEAVE
1436     c save averages
1437     c Phygrow1ave(i,j,k,bi,bj)=Phygrow1ave(i,j,k,bi,bj)+
1438     c & mu1*py1*deltaTclock
1439     c & /float(nsubtime)
1440     c Phygrow2ave(i,j,k,bi,bj)=Phygrow2ave(i,j,k,bi,bj)+
1441     c & mu2*py2*deltaTclock
1442     c & /float(nsubtime)
1443     c Zoograzave(i,j,k,bi,bj)=Zoograzave(i,j,k,bi,bj)+
1444     c & (gampn1*graz1*zo +gampn2*graz2*zo)*
1445     c & deltaTclock/float(nsubtime)
1446     #ifdef GEIDER
1447     Chlave(i,j,k,bi,bj)=Chlave(i,j,k,bi,bj)+
1448     & Chl*dtplankton
1449     #endif
1450     PARave(i,j,k,bi,bj)=PARave(i,j,k,bi,bj)+
1451     & PARl*dtplankton
1452     PPave(i,j,k,bi,bj)=PPave(i,j,k,bi,bj)+
1453     & PP*dtplankton
1454     Nfixave(i,j,k,bi,bj)=Nfixave(i,j,k,bi,bj)+
1455     & Nfix*dtplankton
1456     Denitave(i,j,k,bi,bj)=Denitave(i,j,k,bi,bj)+
1457     & denit*dtplankton
1458     #ifdef WAVES_DIAG_PCHL
1459     do np=1,npmax
1460     Pchlave(i,j,k,bi,bj,np)=Pchlave(i,j,k,bi,bj,np)+
1461     & phychl(np)*dtplankton
1462     enddo
1463     #endif
1464 stephd 1.13 #ifdef DAR_DIAG_PARW
1465     do ilam=1,tlam
1466     PARwave(i,j,k,bi,bj,ilam)=PARwave(i,j,k,bi,bj,ilam)+
1467     & PARw_k(ilam,k)*dtplankton
1468     enddo
1469     do np=1,npmax
1470     chl2cave(i,j,k,bi,bj,np)=chl2cave(i,j,k,bi,bj,np)+
1471     & chl2cl(np)*dtplankton
1472     enddo
1473     #endif
1474 jahn 1.1 #ifdef DAR_DIAG_ACDOM
1475     c print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)
1476     aCDOMave(i,j,k,bi,bj)=aCDOMave(i,j,k,bi,bj)+
1477     & acdom_k(k,darwin_diag_acdom_ilam)*dtplankton
1478     #endif
1479     #ifdef DAR_DIAG_IRR
1480     do ilam = 1,tlam
1481     if (k.EQ.1) then
1482     Edave(i,j,k,bi,bj,ilam)=Edave(i,j,k,bi,bj,ilam)+
1483     & Edwsf(ilam)*dtplankton
1484     Esave(i,j,k,bi,bj,ilam)=Esave(i,j,k,bi,bj,ilam)+
1485     & Eswsf(ilam)*dtplankton
1486     Coj no Eu at surface (yet)
1487     else
1488     Edave(i,j,k,bi,bj,ilam)=Edave(i,j,k,bi,bj,ilam)+
1489     & Edz(ilam,k-1)*dtplankton
1490     Esave(i,j,k,bi,bj,ilam)=Esave(i,j,k,bi,bj,ilam)+
1491     & Esz(ilam,k-1)*dtplankton
1492     Euave(i,j,k,bi,bj,ilam)=Euave(i,j,k,bi,bj,ilam)+
1493     & Euz(ilam,k-1)*dtplankton
1494     endif
1495 jahn 1.10 Estave(i,j,k,bi,bj,ilam)=Estave(i,j,k,bi,bj,ilam)+
1496     & Estop(ilam,k)*dtplankton
1497 jahn 1.1 Eutave(i,j,k,bi,bj,ilam)=Eutave(i,j,k,bi,bj,ilam)+
1498     & Eutop(ilam,k)*dtplankton
1499     enddo
1500     #endif
1501 jahn 1.9 #ifdef DAR_DIAG_IRR_AMPS
1502     do ilam = 1,tlam
1503 jahn 1.10 amp1ave(i,j,k,bi,bj,ilam)=amp1ave(i,j,k,bi,bj,ilam)+
1504     & amp1(ilam,k)*dtplankton
1505     amp2ave(i,j,k,bi,bj,ilam)=amp2ave(i,j,k,bi,bj,ilam)+
1506     & amp2(ilam,k)*dtplankton
1507 jahn 1.9 enddo
1508     #endif
1509 jahn 1.1 #ifdef DAR_DIAG_ABSORP
1510     do ilam = 1,tlam
1511     aave(i,j,k,bi,bj,ilam)=aave(i,j,k,bi,bj,ilam)+
1512     & a_k(k,ilam)*dtplankton
1513     enddo
1514     #endif
1515     #ifdef DAR_DIAG_SCATTER
1516     do ilam = 1,tlam
1517     btave(i,j,k,bi,bj,ilam)=btave(i,j,k,bi,bj,ilam)+
1518     & bt_k(k,ilam)*dtplankton
1519     bbave(i,j,k,bi,bj,ilam)=bbave(i,j,k,bi,bj,ilam)+
1520     & bb_k(k,ilam)*dtplankton
1521     enddo
1522     #endif
1523     #ifdef DAR_DIAG_PART_SCATTER
1524     do ilam = 1,tlam
1525     apartave(i,j,k,bi,bj,ilam)=apartave(i,j,k,bi,bj,ilam)+
1526     & apart_k(k,ilam)*dtplankton
1527     btpartave(i,j,k,bi,bj,ilam)=btpartave(i,j,k,bi,bj,ilam)+
1528     & bpart_k(k,ilam)*dtplankton
1529     bbpartave(i,j,k,bi,bj,ilam)=bbpartave(i,j,k,bi,bj,ilam)+
1530     & bbpart_k(k,ilam)*dtplankton
1531     enddo
1532     #endif
1533 jahn 1.9 #ifdef DAR_RADTRANS
1534     if (k.eq.1) then
1535     rmudave(i,j,bi,bj)=rmudave(i,j,bi,bj)+
1536     & rmud*dtplankton
1537     endif
1538     #endif
1539 stephd 1.13 #ifdef DAR_DIAG_EK
1540     do np=1,npmax
1541     Ekave(i,j,k,bi,bj,np)=Ekave(i,j,k,bi,bj,np)+
1542     & Ekl(np)*dtplankton
1543     EkoverEave(i,j,k,bi,bj,np)=EkoverEave(i,j,k,bi,bj,np)+
1544     & EkoverEl(np)*dtplankton
1545     acclimave(i,j,k,bi,bj,np)=acclimave(i,j,k,bi,bj,np)+
1546     & accliml(np)*dtplankton
1547     do ilam=1,tlam
1548     Ek_nlave(i,j,k,bi,bj,np,ilam)=
1549     & Ek_nlave(i,j,k,bi,bj,np,ilam)+
1550     & Ek_nll(np,ilam)*dtplankton
1551     EkoverE_nlave(i,j,k,bi,bj,np,ilam)=
1552     & EkoverE_nlave(i,j,k,bi,bj,np,ilam)+
1553     & EkoverE_nll(np,ilam)*dtplankton
1554     enddo
1555     enddo
1556     #endif
1557 jahn 1.1 #ifdef DAR_DIAG_RSTAR
1558     do np=1,npmax
1559     Rstarave(i,j,k,bi,bj,np)=Rstarave(i,j,k,bi,bj,np)+
1560     & Rstarl(np)*dtplankton
1561     RNstarave(i,j,k,bi,bj,np)=RNstarave(i,j,k,bi,bj,np)+
1562     & RNstarl(np)*dtplankton
1563     enddo
1564     #endif
1565     #ifdef DAR_DIAG_DIVER
1566     Diver1ave(i,j,k,bi,bj)=Diver1ave(i,j,k,bi,bj)+
1567     & Diver1(i,j,k)*dtplankton
1568     Diver2ave(i,j,k,bi,bj)=Diver2ave(i,j,k,bi,bj)+
1569     & Diver2(i,j,k)*dtplankton
1570     Diver3ave(i,j,k,bi,bj)=Diver3ave(i,j,k,bi,bj)+
1571     & Diver3(i,j,k)*dtplankton
1572     Diver4ave(i,j,k,bi,bj)=Diver4ave(i,j,k,bi,bj)+
1573     & Diver4(i,j,k)*dtplankton
1574     #endif
1575     #ifdef DAR_DIAG_GROW
1576     do np=1,npmax
1577     Growave(i,j,k,bi,bj,np)=Growave(i,j,k,bi,bj,np)+
1578     & Growl(np)*dtplankton
1579     Growsqave(i,j,k,bi,bj,np)=Growsqave(i,j,k,bi,bj,np)+
1580     & Growsql(np)*dtplankton
1581     enddo
1582     #endif
1583    
1584     #ifdef ALLOW_DIAZ
1585     #ifdef DAR_DIAG_NFIXP
1586     do np=1,npmax
1587     NfixPave(i,j,k,bi,bj,np)=NfixPave(i,j,k,bi,bj,np)+
1588     & NfixPl(np)*dtplankton
1589     enddo
1590     #endif
1591     #endif
1592     #endif
1593    
1594     #ifdef ALLOW_CARBON
1595     if (k.eq.1) then
1596     SURave(i,j,bi,bj) =SURave(i,j,bi,bj)+
1597 jahn 1.2 & flxCO2(i,j)*dtplankton
1598 jahn 1.1 SURCave(i,j,bi,bj) =SURCave(i,j,bi,bj)+
1599     & FluxCO2(i,j,bi,bj)*dtplankton
1600     SUROave(i,j,bi,bj) =SUROave(i,j,bi,bj)+
1601 jahn 1.2 & flxO2(i,j)*dtplankton
1602 stephd 1.15 endif
1603     #ifdef pH_3D
1604     pCO2ave(i,j,k,bi,bj) =pCO2ave(i,j,k,bi,bj)+
1605     & pCO2(i,j,k,bi,bj)*dtplankton
1606     pHave(i,j,k,bi,bj) =pHave(i,j,k,bi,bj)+
1607     & pH(i,j,k,bi,bj)*dtplankton
1608     #else
1609     if (k.eq.1) then
1610 jahn 1.1 pCO2ave(i,j,bi,bj) =pCO2ave(i,j,bi,bj)+
1611     & pCO2(i,j,bi,bj)*dtplankton
1612     pHave(i,j,bi,bj) =pHave(i,j,bi,bj)+
1613     & pH(i,j,bi,bj)*dtplankton
1614     endif
1615     #endif
1616 stephd 1.15 #endif
1617 jahn 1.1 endif
1618     c end if hFac>0
1619    
1620     enddo ! k
1621     c end layer loop
1622     c
1623    
1624     ENDDO ! i
1625     ENDDO ! j
1626    
1627     #ifdef ALLOW_PAR_DAY
1628     C 1 <-> 2
1629     PARiaccum = 3 - PARiprev
1630    
1631     DO k=1,nR
1632     DO j=1,sNy
1633     DO i=1,sNx
1634     PARday(i,j,k,bi,bj,PARiaccum) =
1635     & PARday(i,j,k,bi,bj,PARiaccum) + PAR(i,j,k)
1636     ENDDO
1637     ENDDO
1638     ENDDO
1639    
1640     phase = 0. _d 0
1641     itistime = DIFF_PHASE_MULTIPLE( phase, darwin_PARavPeriod,
1642     & newtime, dtsubtime)
1643    
1644     IF ( itistime ) THEN
1645     C compute average
1646     nav = darwin_PARnav
1647     IF (newtime - baseTime .LT. darwin_PARavPeriod) THEN
1648     C incomplete period at beginning of run
1649     nav = NINT((newtime-baseTime)/dtsubtime)
1650     ENDIF
1651     DO k=1,nR
1652     DO j=1,sNy
1653     DO i=1,sNx
1654     PARday(i,j,k,bi,bj,PARiaccum) =
1655     & PARday(i,j,k,bi,bj,PARiaccum) / nav
1656     ENDDO
1657     ENDDO
1658     ENDDO
1659     C reset the other slot for averaging
1660     DO k=1,nR
1661     DO j=1,sNy
1662     DO i=1,sNx
1663     PARday(i,j,k,bi,bj,PARiprev) = 0. _d 0
1664     ENDDO
1665     ENDDO
1666     ENDDO
1667     ENDIF
1668     C itistime
1669     #endif
1670    
1671 jahn 1.11 #ifdef DAR_CHECK_IRR_CONT
1672 jahn 1.12 i = myXGlobalLo-1+(bi-1)*sNx+idiscEs
1673     j = myYGlobalLo-1+(bj-1)*sNy+jdiscEs
1674     write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Es disc',
1675     & i,j,kdiscEs,ldiscEs,discEs
1676     i = myXGlobalLo-1+(bi-1)*sNx+idiscEu
1677     j = myYGlobalLo-1+(bj-1)*sNy+jdiscEu
1678     write(6,'(I4.4,X,A,4(X,I4),1PE24.16)')myProcId,'max Eu disc',
1679     & i,j,kdiscEu,ldiscEu,discEu
1680 jahn 1.11 #endif
1681    
1682 jahn 1.1 COJ fill diagnostics
1683     #ifdef ALLOW_DIAGNOSTICS
1684     IF ( useDiagnostics ) THEN
1685     diagname = ' '
1686     WRITE(diagname,'(A8)') 'PAR '
1687     CALL DIAGNOSTICS_FILL( PAR(1-Olx,1-Oly,1), diagname,
1688     & 0,Nr,2,bi,bj,myThid )
1689     WRITE(diagname,'(A8)') 'PP '
1690     CALL DIAGNOSTICS_FILL( PParr(1-Olx,1-Oly,1), diagname,
1691     & 0,Nr,2,bi,bj,myThid )
1692     WRITE(diagname,'(A8)') 'Nfix '
1693     CALL DIAGNOSTICS_FILL( Nfixarr(1-Olx,1-Oly,1), diagname,
1694     & 0,Nr,2,bi,bj,myThid )
1695     c ANNA_TAVE
1696     #ifdef WAVES_DIAG_PCHL
1697     DO np=1,MIN(99,npmax)
1698     WRITE(diagname,'(A5,I2.2,A1)') 'Pchl',np,' '
1699     CALL DIAGNOSTICS_FILL( Pchlarr(1-Olx,1-Oly,1,np), diagname,
1700     & 0,Nr,2,bi,bj,myThid )
1701     ENDDO
1702     #endif
1703     c ANNA end TAVE
1704     #ifdef DAR_DIAG_RSTAR
1705     DO np=1,MIN(99,npmax)
1706     WRITE(diagname,'(A5,I2.2,A1)') 'Rstar',np,' '
1707     CALL DIAGNOSTICS_FILL( Rstararr(1-Olx,1-Oly,1,np), diagname,
1708     & 0,Nr,2,bi,bj,myThid )
1709     ENDDO
1710     #endif
1711     #ifdef DAR_DIAG_DIVER
1712     WRITE(diagname,'(A8)') 'Diver1 '
1713     CALL DIAGNOSTICS_FILL( Diver1(1-Olx,1-Oly,1), diagname,
1714     & 0,Nr,2,bi,bj,myThid )
1715     WRITE(diagname,'(A8)') 'Diver2 '
1716     CALL DIAGNOSTICS_FILL( Diver2(1-Olx,1-Oly,1), diagname,
1717     & 0,Nr,2,bi,bj,myThid )
1718     WRITE(diagname,'(A8)') 'Diver3 '
1719     CALL DIAGNOSTICS_FILL( Diver3(1-Olx,1-Oly,1), diagname,
1720     & 0,Nr,2,bi,bj,myThid )
1721     WRITE(diagname,'(A8)') 'Diver4 '
1722     CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
1723     & 0,Nr,2,bi,bj,myThid )
1724 jahn 1.5 WRITE(diagname,'(A8)') 'Shannon '
1725     CALL DIAGNOSTICS_FILL( Shannon(1-Olx,1-Oly,1), diagname,
1726     & 0,Nr,2,bi,bj,myThid )
1727     WRITE(diagname,'(A8)') 'Simpson '
1728     CALL DIAGNOSTICS_FILL( Simpson(1-Olx,1-Oly,1), diagname,
1729     & 0,Nr,2,bi,bj,myThid )
1730 jahn 1.1 #endif
1731     #ifdef ALLOW_DIAZ
1732     #ifdef DAR_DIAG_NFIXP
1733     DO np=1,MIN(99,npmax)
1734     WRITE(diagname,'(A5,I2.2,A1)') 'NfixP',np,' '
1735     CALL DIAGNOSTICS_FILL( NfixParr(1-Olx,1-Oly,1,np), diagname,
1736     & 0,Nr,2,bi,bj,myThid )
1737     ENDDO
1738     #endif
1739     #endif
1740     #ifdef DAR_DIAG_CHL
1741     CALL DIAGNOSTICS_FILL( GeiderChlarr(1-Olx,1-Oly,1), 'ChlGeide',
1742     & 0,Nr,2,bi,bj,myThid )
1743     CALL DIAGNOSTICS_FILL( GeiderChl2Carr(1-Olx,1-Oly,1),'Chl2CGei',
1744     & 0,Nr,2,bi,bj,myThid )
1745     CALL DIAGNOSTICS_FILL( DoneyChlarr(1-Olx,1-Oly,1), 'ChlDoney',
1746     & 0,Nr,2,bi,bj,myThid )
1747     CALL DIAGNOSTICS_FILL( DoneyChl2Carr(1-Olx,1-Oly,1), 'Chl2CDon',
1748     & 0,Nr,2,bi,bj,myThid )
1749     CALL DIAGNOSTICS_FILL( CloernChlarr(1-Olx,1-Oly,1), 'ChlCloer',
1750     & 0,Nr,2,bi,bj,myThid )
1751     CALL DIAGNOSTICS_FILL( CloernChl2Carr(1-Olx,1-Oly,1),'Chl2CClo',
1752     & 0,Nr,2,bi,bj,myThid )
1753     #endif
1754     #ifdef ALLOW_CARBON
1755 jahn 1.2 CALL DIAGNOSTICS_FILL( flxCO2(1-Olx,1-Oly), 'DICTFLX ',
1756 jahn 1.1 & 0,1,2,bi,bj,myThid )
1757     CALL DIAGNOSTICS_FILL( FluxCO2(1-Olx,1-Oly,bi,bj), 'DICCFLX ',
1758     & 0,1,2,bi,bj,myThid )
1759 jahn 1.2 CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly), 'DICOFLX ',
1760 jahn 1.1 & 0,1,2,bi,bj,myThid )
1761 stephd 1.15 #ifdef pH_3D
1762     CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,1,bi,bj), 'DICPCO2 ',
1763     & 0,Nr,2,bi,bj,myThid )
1764     CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,1,bi,bj), 'DICPHAV ',
1765     & 0,Nr,2,bi,bj,myThid )
1766     #else
1767 jahn 1.1 CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',
1768     & 0,1,2,bi,bj,myThid )
1769     CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,bi,bj), 'DICPHAV ',
1770     & 0,1,2,bi,bj,myThid )
1771 stephd 1.15 #endif
1772 jahn 1.1 #endif /* ALLOW_CARBON */
1773     ENDIF
1774     #endif /* ALLOW_DIAGNOSTICS */
1775     COJ
1776    
1777     c determine iron partitioning - solve for free iron
1778     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
1779     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
1780     & myIter, mythid)
1781     c
1782     #ifdef ALLOW_TIMEAVE
1783     c save averages
1784 jahn 1.16 dar_timeave(bi,bj) = dar_timeave(bi,bj) + dtplankton
1785 jahn 1.1 #ifdef ALLOW_CARBON
1786 jahn 1.16 dic_timeave(bi,bj) = dic_timeave(bi,bj) + dtplankton
1787 jahn 1.1 #endif
1788     #endif
1789     c
1790     c -----------------------------------------------------
1791     ENDDO ! it
1792     c -----------------------------------------------------
1793     c end of bio-chemical time loop
1794     c
1795     RETURN
1796     END
1797     #endif /*MONOD*/
1798     #endif /*ALLOW_PTRACERS*/
1799    
1800     C============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22