/[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.6 - (hide annotations) (download)
Thu May 31 21:08:25 2012 UTC (13 years, 2 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt63n_20120604
Changes since 1.5: +31 -1 lines
o add CDOM-like tracer (#define ALLOW_CDOM)

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

  ViewVC Help
Powered by ViewVC 1.1.22