/[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.1 - (hide annotations) (download)
Wed Apr 13 18:56:25 2011 UTC (14 years, 3 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt62v_20110413, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt62w_20110426, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_baseline, ctrb_darwin2_ckpt62z_20110622
darwin2 initial checkin

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

  ViewVC Help
Powered by ViewVC 1.1.22