/[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.18 - (hide annotations) (download)
Fri Jul 24 21:47:26 2015 UTC (10 years ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt65r_20151221
Changes since 1.17: +4 -2 lines
make it work with useDiagnostics false

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

  ViewVC Help
Powered by ViewVC 1.1.22