/[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.7 - (hide annotations) (download)
Fri Jun 29 20:41:59 2012 UTC (13 years ago) by stephd
Branch: MAIN
Changes since 1.6: +18 -1 lines
o add potential for CO2 limitation to growth rate -- commented out for now
  (need to check monod_forcing and monod_plankton to make it work)

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

  ViewVC Help
Powered by ViewVC 1.1.22