/[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.20 - (hide annotations) (download)
Fri Apr 1 00:02:17 2016 UTC (9 years, 3 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt66m_20171213, HEAD
Changes since 1.19: +3 -3 lines
reduce light due to ice also with radtrans
since OASIM fields do not include its effect

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

  ViewVC Help
Powered by ViewVC 1.1.22