/[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.4 - (hide annotations) (download)
Mon Nov 7 21:55:11 2011 UTC (13 years, 8 months ago) by jahn
Branch: MAIN
Changes since 1.3: +9 -1 lines
initialize diagnostic chlorophyll arrays

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

  ViewVC Help
Powered by ViewVC 1.1.22