/[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.5 - (hide annotations) (download)
Wed Nov 9 23:34:17 2011 UTC (13 years, 8 months ago) by jahn
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63k_20120317
Changes since 1.4: +27 -1 lines
add Shannon and Simpson diversity indices

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

  ViewVC Help
Powered by ViewVC 1.1.22