/[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.8 - (hide annotations) (download)
Sat Jun 30 19:23:04 2012 UTC (13 years ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63p_20120707
Changes since 1.7: +4 -7 lines
avoid double inclusion

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

  ViewVC Help
Powered by ViewVC 1.1.22