/[MITgcm]/MITgcm_contrib/dcarroll/highres_darwin/code/darwin_forcing.F
ViewVC logotype

Annotation of /MITgcm_contrib/dcarroll/highres_darwin/code/darwin_forcing.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Sun Sep 22 21:23:46 2019 UTC (5 years, 10 months ago) by dcarroll
Branch: MAIN
CVS Tags: HEAD
Initial check in of high resolution Darwin simulation code

1 dcarroll 1.1 C $Header: /u/gcmpack/MITgcm_contrib/ecco_darwin/v4_llc270/code_darwin/darwin_forcing.F,v 1.17 2019/09/16 15:25:49 dcarroll Exp $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5     #include "PTRACERS_OPTIONS.h"
6     #include "DARWIN_OPTIONS.h"
7    
8     #ifdef ALLOW_PTRACERS
9     #ifdef ALLOW_DARWIN
10    
11     c=============================================================
12     c subroutine DARWIN_forcing
13     c step forward bio-chemical tracers in time
14     C==============================================================
15     SUBROUTINE DARWIN_Forcing(
16     U Ptr,
17     I bi,bj,imin,imax,jmin,jmax,
18     I myIter,myTime,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 "DARWIN_SIZE.h"
34     #include "DARWIN.h"
35     #include "DARWIN_IO.h"
36     #include "DARWIN_FLUX.h"
37     #include "DARWIN_FIELDS.h"
38     #include "GGL90.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    
48     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     c iron partitioning
74     _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
75     c some working variables
76     _RL sumpy
77     _RL sumpyup
78     c light variables
79     _RL PAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
80     _RL sfac(1-OLy:sNy+OLy)
81     _RL atten,lite
82     _RL newtime ! for sub-timestepping
83     _RL runtim ! time from tracer initialization
84    
85    
86     c ANNA define variables for wavebands
87     #ifdef WAVEBANDS
88     integer ilam
89     _RL PARw_k(tlam,Nr)
90     _RL PARwup(tlam)
91     _RL acdom_k(Nr,tlam)
92     #ifdef DAR_RADTRANS
93     integer iday,iyr,imon,isec,lp,wd,mydate(4)
94     _RL Edwsf(tlam),Eswsf(tlam)
95     _RL Edz(tlam,Nr),Esz(tlam,Nr),Euz(tlam,Nr),Eutop(tlam,Nr)
96     _RL tirrq(nr)
97     _RL tirrwq(tlam,nr)
98     _RL solz
99     _RL rmud
100     _RL actot,bctot,bbctot
101     _RL apart_k(Nr,tlam),bpart_k(Nr,tlam),bbpart_k(Nr,tlam)
102     _RL bt_k(Nr,tlam), bb_k(Nr,tlam)
103     #else
104     _RL PARwdn(tlam)
105     #endif
106     C always need for diagnostics
107     _RL a_k(Nr,tlam)
108     #endif /* WAVEBANDS */
109    
110    
111     #ifdef DAR_DIAG_DIVER
112     _RL Diver1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
113     _RL Diver2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114     _RL Diver3(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115     _RL Diver4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116    
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     _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    
257     _RL calcium(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
258     _RL KspTP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
259     #ifdef ADKINS_SURF_FLUX
260     _RL dcal
261     #endif
262     #ifdef ALLOW_SED_DISS_FLUX
263     c sediment-to-ocean fluxes
264     _RL DICSedFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
265     _RL ALKSedFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
266     _RL BBLDiffusionCoeffLoc
267     _RL BBLThicknessLoc
268     _RL RFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
269     _RL CO3Sw(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
270     _RL CO3Sed(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
271     #endif /* ALLOW_SED_DISS_FLUX */
272    
273     _RL t
274     _RL s
275     _RL ta
276     _RL pt
277     _RL sit
278     _RL tk
279     _RL tk100
280     _RL tk1002
281     _RL dlogtk
282     _RL sqrtis
283     _RL sqrts
284     _RL s15
285     _RL scl
286     _RL x1
287     _RL x2
288     _RL s2
289     _RL xacc
290     _RL invtk
291     _RL is
292     _RL is2
293     _RL bdepth
294     _RL cdepth
295     _RL pressc
296     _RL Ksp_T_Calc
297     _RL xvalue
298     _RL zdum
299     _RL tmpa1
300     _RL tmpa2
301     _RL tmpa3
302     _RL logKspc
303     _RL dv
304     _RL dk
305     _RL pfactor
306     _RL bigR
307    
308     _RL pCO2SolverTemp
309     _RL pCO2SolverSal
310     _RL pCO2SolverDic
311     _RL pCO2SolverPo4
312     _RL pCO2SolverSi
313     _RL pCO2SolverAlk
314    
315     #ifdef CO2_FLUX_BUDGET
316     _RL pHBudget1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
317     _RL pCO2Budget1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
318     _RL CO3Budget1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
319    
320     _RL baselinePH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
321     _RL baselinePCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
322     _RL baselineCO3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
323    
324     _RL pCO2Theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
325     _RL pCO2Salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
326     _RL pCO2Alk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
327     _RL pCO2Dic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
328    
329     _RL deltaTheta(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
330     _RL deltaSalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
331     _RL deltaAlk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
332     _RL deltaDic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
333     _RL deltaApCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
334    
335     _RL deltaDic_theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
336     _RL deltaDic_salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
337     _RL deltaDic_alk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
338     _RL deltaDic_apCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
339     _RL deltaDic_bio(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
340     _RL deltaDic_residual(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
341    
342     _RL mixingDepthKLev(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
343     _RL mixingDepth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
344    
345     _RL dCO2Flux_theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
346     _RL dCO2Flux_salt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
347     _RL dCO2Flux_alk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
348     _RL dCO2Flux_dic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
349     _RL dCO2Flux_apCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
350     _RL dCO2Flux_residual(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
351     _RL dCO2Flux_bio(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
352     _RL dCO2Flux_circ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
353     #endif /* CO2_FLUX_BUDGET */
354     #endif
355    
356     _RL tot_Nfix
357     _RL tmp
358     _RL phytmp, chltmp
359    
360     INTEGER i,j,k,it, ktmp
361     INTEGER np, nz, np2, npsave
362     INTEGER debug
363     CHARACTER*8 diagname
364    
365     DO j=1-OLy,sNy+OLy
366     DO i=1-OLx,sNx+OLx
367     do k=1,Nr
368     freefe(i,j,k)=0. _d 0
369     PAR(i,j,k) = 0. _d 0
370     #ifdef DAR_DIAG_DIVER
371     Diver1(i,j,k)=0. _d 0
372     Diver2(i,j,k)=0. _d 0
373     Diver3(i,j,k)=0. _d 0
374     Diver4(i,j,k)=0. _d 0
375     #endif
376    
377     #ifdef ALLOW_DIAGNOSTICS
378     COJ for diagnostics
379     PParr(i,j,k) = 0. _d 0
380     Nfixarr(i,j,k) = 0. _d 0
381     #ifdef DAR_DIAG_CHL
382     GeiderChlarr(i,j,k) = 0. _d 0
383     GeiderChl2Carr(i,j,k) = 0. _d 0
384     DoneyChlarr(i,j,k) = 0. _d 0
385     DoneyChl2Carr(i,j,k) = 0. _d 0
386     CloernChlarr(i,j,k) = 0. _d 0
387     CloernChl2Carr(i,j,k) = 0. _d 0
388     #endif
389     c ANNA_TAVE
390     #ifdef WAVES_DIAG_PCHL
391     DO np=1,npmax
392     Pchlarr(i,j,k,np) = 0. _d 0
393     ENDDO
394     #endif
395     c ANNA end TAVE
396     #ifdef DAR_DIAG_RSTAR
397     DO np=1,npmax
398     Rstararr(i,j,k,np) = 0. _d 0
399     ENDDO
400     #endif
401     COJ
402     #ifdef ALLOW_DIAZ
403     #ifdef DAR_DIAG_NFIXP
404     DO np=1,npmax
405     NfixParr(i,j,k,np) = 0. _d 0
406     ENDDO
407     #endif
408     #endif
409     #endif
410     enddo
411     ENDDO
412     ENDDO
413     c
414     c bio-chemical time loop
415     c--------------------------------------------------
416     DO it=1,nsubtime
417     c -------------------------------------------------
418     tot_Nfix=0. _d 0
419     COJ cannot use dfloat because of adjoint
420     COJ division will be double precision anyway because of dTtracerLev
421     newtime=myTime-dTtracerLev(1)+
422     & float(it)*dTtracerLev(1)/float(nsubtime)
423     c print*,'it ',it,newtime,nsubtime,myTime
424     runtim=myTime-float(PTRACERS_Iter0)*dTtracerLev(1)
425    
426     c determine iron partitioning - solve for free iron
427     c ---------------------------
428     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
429     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
430     & myIter, mythid)
431     c --------------------------
432     #ifdef ALLOW_CARBON
433    
434     #ifdef CO2_FLUX_BUDGET
435     C store values from previous timestep
436     DO j=jmin,jmax
437     DO i=imin,imax
438     pHBudget1(i,j) = pH(i,j,bi,bj)
439     pCO2Budget1(i,j) = pCO2(i,j,bi,bj)
440     CO3Budget1(i,j) = CO3(i,j,bi,bj)
441     ENDDO
442     ENDDO
443     #endif /* CO2_FLUX_BUDGET */
444    
445     C compute baseline pCO2 and air-sea CO2 flux
446     call dic_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
447     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
448     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
449     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
450     & flxCO2,
451     & bi,bj,imin,imax,jmin,jmax,
452     & myIter,myTime,myThid)
453    
454     #ifdef CO2_FLUX_BUDGET
455     C store values from current timestep as non-perturbed baseline
456     DO j=jmin,jmax
457     DO i=imin,imax
458     baselinePH(i,j) = pH(i,j,bi,bj)
459     baselinePCO2(i,j) = pCO2(i,j,bi,bj)
460     baselineCO3(i,j) = CO3(i,j,bi,bj)
461     ENDDO
462     ENDDO
463    
464     C set pH to value from previous timestep
465     DO j=jmin,jmax
466     DO i=imin,imax
467     pH(i,j,bi,bj) = pHBudget1(i,j)
468     ENDDO
469     ENDDO
470    
471     C deltaPCO2 due to theta pertubation
472     call dic_budgetTheta(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
473     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
474     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
475     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
476     & deltaTheta,bi,bj,imin,imax,jmin,jmax,
477     & myIter,myTime,myThid)
478    
479     C set pH to value from previous timestep and store pCO2
480     DO j=jmin,jmax
481     DO i=imin,imax
482     pH(i,j,bi,bj) = pHBudget1(i,j)
483     pCO2Theta(i,j) = pCO2(i,j,bi,bj)
484     ENDDO
485     ENDDO
486    
487     C deltaPCO2 due to salinity perturbation
488     call dic_budgetSalt(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
489     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
490     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
491     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
492     & deltaSalt,bi,bj,imin,imax,jmin,jmax,
493     & myIter,myTime,myThid)
494    
495     C set pH to value from previous timestep and store pCO2
496     DO j=jmin,jmax
497     DO i=imin,imax
498     pH(i,j,bi,bj) = pHBudget1(i,j)
499     pCO2Salt(i,j) = pCO2(i,j,bi,bj)
500     ENDDO
501     ENDDO
502    
503     C deltaPCO2 due to alkalinity perturbation
504     call dic_budgetAlk(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
505     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
506     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
507     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
508     & deltaAlk,bi,bj,imin,imax,jmin,jmax,
509     & myIter,myTime,myThid)
510    
511     C set pH to value from previous timestep and store pCO2
512     DO j=jmin,jmax
513     DO i=imin,imax
514     pH(i,j,bi,bj) = pHBudget1(i,j)
515     pCO2Alk(i,j) = pCO2(i,j,bi,bj)
516     ENDDO
517     ENDDO
518    
519     C deltaPCO2 due to DIC perturbation
520     call dic_budgetDic(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
521     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
522     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
523     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
524     & deltaDic,bi,bj,imin,imax,jmin,jmax,
525     & myIter,myTime,myThid)
526    
527     C set pH to value from previous timestep and store pCO2
528     DO j=jmin,jmax
529     DO i=imin,imax
530     pH(i,j,bi,bj) = pHBudget1(i,j)
531     pCO2Dic(i,j) = pCO2(i,j,bi,bj)
532     ENDDO
533     ENDDO
534    
535     C compute deltaApCO2
536     call dic_budgetApCO2(Ptr(1-OLx,1-OLy,1,bi,bj,iDIC),
537     & Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
538     & Ptr(1-OLx,1-OLy,1,bi,bj,iPO4),
539     & Ptr(1-OLx,1-OLy,1,bi,bj,iSi),
540     & deltaApCO2,bi,bj,imin,imax,jmin,jmax,
541     & myIter,myTime,myThid)
542    
543     C reset to baseline values
544     DO j=jmin,jmax
545     DO i=imin,imax
546     pH(i,j,bi,bj) = baselinePH(i,j)
547     pCO2(i,j,bi,bj) = baselinePCO2(i,j)
548     CO3(i,j,bi,bj) = baselineCO3(i,j)
549     ENDDO
550     ENDDO
551    
552     DO j=jmin,jmax
553     DO i=imin,imax
554    
555     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
556    
557     C compute delta DIC terms for each budget component (mmol C m^-3)
558     deltaDic_theta(i,j) =
559     & (((pCO2Theta(i,j) - baselinePCO2(i,j)) /
560     & budgetPert) /
561     & ((pCO2Dic(i,j) - baselinePCO2(i,j)) /
562     & budgetPert)) *
563     & deltaTheta(i,j)
564    
565     deltaDic_salt(i,j) =
566     & (((pCO2Salt(i,j) - baselinePCO2(i,j)) /
567     & budgetPert) /
568     & ((pCO2Dic(i,j) - baselinePCO2(i,j)) /
569     & budgetPert)) *
570     & deltaSalt(i,j)
571    
572     deltaDic_alk(i,j) =
573     & (((pCO2Alk(i,j) - baselinePCO2(i,j)) /
574     & budgetPert) /
575     & ((pCO2Dic(i,j) - baselinePCO2(i,j)) /
576     & budgetPert)) *
577     & deltaAlk(i,j)
578    
579     deltaDic_apCO2(i,j) =
580     & (budgetPert / (pCO2Dic(i,j) - baselinePCO2(i,j))) *
581     & deltaApCO2(i,j)
582    
583     deltaDic_residual(i,j) = deltaDic(i,j) -
584     & (deltaDic_theta(i,j) + deltaDic_salt(i,j) +
585     & deltaDic_alk(i,j) + deltaDic_apCO2(i,j))
586    
587     else
588    
589     deltaDic_theta(i,j) = 0. _d 0
590     deltaDic_salt(i,j) = 0. _d 0
591     deltaDic_alk(i,j) = 0. _d 0
592     deltaDic_apCO2(i,j) = 0. _d 0
593     deltaDic_residual(i,j) = 0. _d 0
594    
595     endif
596    
597     ENDDO
598     ENDDO
599     #endif /* CO2_FLUX_BUDGET */
600     c air-sea flux of O2
601     call dic_o2_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iO2),
602     & flxO2,
603     & bi,bj,imin,imax,jmin,jmax,
604     & myIter,myTime,myThid)
605     c dilution of alkalinity
606     call dic_alk_surfforcing(Ptr(1-OLx,1-OLy,1,bi,bj,iALK),
607     & flxALK,
608     & bi,bj,imin,imax,jmin,jmax,
609     & myIter,myTime,myThid)
610     #endif
611    
612    
613     c find light in each grid cell
614     c ---------------------------
615     c determine incident light
616     #ifndef READ_PAR
617     #ifndef USE_QSW
618     DO j=1-OLy,sNy+OLy
619     sfac(j)=0. _d 0
620     ENDDO
621     call darwin_insol(newTime,sfac,bj)
622     #endif /* not USE_QSW */
623     #endif /* not READ_PAR */
624    
625     #ifdef ALLOW_PAR_DAY
626     C find out which slot of PARday has previous day's average
627     dtsubtime = dTtracerLev(1)/float(nsubtime)
628     C running index of averaging period
629     C myTime has already been incremented in this iteration,
630     C go back half a substep to avoid roundoff problems
631     iperiod = FLOOR((newtime-0.5 _d 0*dtsubtime)
632     & /darwin_PARavPeriod)
633     C 0 -> 1, 1->2, 2->0, ...
634     PARiprev = MOD(iperiod, 2) + 1
635    
636     #ifdef ALLOW_DIAGNOSTICS
637     C always fill; this will be the same during PARavPeriod, but this
638     C way it won't blow up for weird diagnostics periods.
639     C we fill before updating, so the diag is the one used in this time
640     C step
641     IF ( useDiagnostics ) THEN
642     CALL DIAGNOSTICS_FILL(
643     & PARday(1-Olx,1-Oly,1,bi,bj,PARiprev),'PARday ',
644     & 0,Nr,2,bi,bj,myThid )
645     ENDIF
646     #endif
647     #endif /* ALLOW_PAR_DAY */
648    
649     #ifdef DAR_RADTRANS
650     #ifndef DAR_RADTRANS_USE_MODEL_CALENDAR
651     #ifdef ALLOW_CAL
652     C get current date and time of day: iyr/imon/iday+isec
653     CALL CAL_GETDATE( myIter, newtime, mydate, mythid )
654     CALL CAL_CONVDATE( mydate,iyr,imon,iday,isec,lp,wd,mythid )
655     #else
656     STOP 'need cal package or DAR_RADTRANS_USE_MODEL_CALENDAR'
657     #endif
658     #endif
659     #endif
660    
661     C.................................................................
662     C.................................................................
663    
664    
665     C ========================== i,j loops =================================
666     DO j=1,sNy
667     DO i=1,sNx
668    
669     c ------------ these are convenient ------------------------------------
670     DO k=1,Nr
671     part_k(k) = max(Ptr(i,j,k,bi,bj,iPOP),0. _d 0)
672     DO np = 1,npmax
673     Phy_k(np,k) = max(Ptr(i,j,k,bi,bj,iPhy+np-1),0. _d 0)
674     #ifdef GEIDER
675     #ifdef DYNAMIC_CHL
676     phychl_k(np,k) = max(Ptr(i,j,k,bi,bj,iChl+np-1),0. _d 0)
677     #else
678     phychl_k(np,k) = max(Chl_phy(i,j,k,bi,bj,np), 0. _d 0)
679     #endif
680     #endif
681     ENDDO
682     ENDDO
683    
684     c ------------ GET CDOM_k FOR WAVEBANDS_3D and RADTRANS ----------------
685     #ifdef WAVEBANDS
686     #if defined(DAR_CALC_ACDOM) || defined(DAR_RADTRANS)
687     call darwin_acdom(phychl_k,aphy_chl,aw,
688     O acdom_k,
689     I myThid)
690     #else
691     DO k=1,Nr
692     DO ilam = 1,tlam
693     acdom_k(k,ilam) = acdom(ilam)
694     ENDDO
695     ENDDO
696     #endif /* DAR_CALC_ACDOM or DAR_RADTRANS */
697     #endif /* WAVEBANDS */
698    
699     c ------------ GET INCIDENT NON-SPECTRAL LIGHT -------------------------
700     #if !(defined(WAVEBANDS) && defined(OASIM))
701     #ifdef READ_PAR
702    
703     lite = sur_par(i,j,bi,bj)
704    
705     #else /* not READ_PAR */
706     #ifdef USE_QSW
707    
708     #ifdef ALLOW_LONGSTEP
709     Qswlocal=LS_Qsw(i,j,bi,bj)
710     #else
711     Qswlocal=Qsw(i,j,bi,bj)
712     #endif
713     lite = -parfrac*Qswlocal*parconv*maskC(i,j,1,bi,bj)
714    
715     #else /* not USE_QSW */
716    
717     C convert W/m2 to uEin/s/m2
718     lite = sfac(j)*parconv*maskC(i,j,1,bi,bj)
719    
720     #endif /* not USE_QSW */
721     #endif /* not READ_PAR */
722    
723     c take ice coverage into account
724     c unless already done in seaice package
725     #if !(defined (ALLOW_SEAICE) && defined (USE_QSW))
726     lite = lite*(1. _d 0-fice(i,j,bi,bj))
727     #endif
728     #endif /* not(WAVEBANDS and OASIM) */
729    
730     c ------------ LIGHT ATTENUATION: --------------------------------------
731     #ifndef WAVEBANDS
732     c ------------ SINGLE-BAND ATTENUATION ---------------------------------
733     atten=0. _d 0
734     do k=1,Nr
735     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
736     sumpyup = sumpy
737     sumpy = 0. _d 0
738     do np=1,npmax
739     #ifdef GEIDER
740     sumpy = sumpy + phychl_k(np,k)
741     #else
742     sumpy = sumpy + Phy_k(np,k)
743     #endif
744     enddo
745     atten= atten + (k0 + kc*sumpy)*5. _d -1*drF(k)
746     if (k.gt.1)then
747     atten = atten + (k0+kc*sumpyup)*5. _d -1*drF(k-1)
748     endif
749     PAR(i,j,k) = lite*exp(-atten)
750     endif
751     enddo
752    
753     #else /* WAVEBANDS */
754     #ifndef DAR_RADTRANS
755     c ------------ WAVEBANDS W/O RADTRANS ----------------------------------
756     do ilam = 1,tlam
757     #ifdef OASIM
758     c add direct and diffuse, convert to uEin/m2/s/nm
759     PARwup(ilam) = WtouEins(ilam)*(oasim_ed(i,j,ilam,bi,bj)+
760     & oasim_es(i,j,ilam,bi,bj))
761     c and take ice fraction into account
762     c PARwup(ilam) = PARwup(ilam)*(1 _d 0 - fice(i,j,bi,bj))
763     #else
764     c sf is per nm; convert to per waveband
765     PARwup(ilam) = wb_width(ilam)*sf(ilam)*lite
766     #endif
767     enddo
768    
769     do k=1,Nr
770     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
771     do ilam = 1,tlam
772     sumpy = 0.
773     do np = 1,npmax
774     c get total attenuation (absorption) by phyto at each wavelength
775     sumpy = sumpy + (phychl_k(np,k)*aphy_chl(np,ilam))
776     enddo
777     c for diagnostic
778     a_k(k,ilam) = aw(ilam) + sumpy + acdom_k(k,ilam)
779     atten = a_k(k,ilam)*drF(k)
780     PARwdn(ilam) = PARwup(ilam)*exp(-atten)
781     enddo
782    
783     c find for the midpoint of the gridcell (gridcell mean)
784     do ilam = 1,tlam
785     C PARw_k(ilam,k)=exp((log(PARwup(ilam))+log(PARwdn(ilam)))*0.5)
786     PARw_k(ilam,k)=sqrt(PARwup(ilam)*PARwdn(ilam))
787     enddo
788    
789     c cycle
790     do ilam=1,tlam
791     PARwup(ilam) = PARwdn(ilam)
792     enddo
793     else
794     do ilam=1,tlam
795     PARw_k(ilam,k) = 0. _d 0
796     enddo
797     endif
798    
799     c sum wavebands for total PAR at the mid point of the gridcell (PARl)
800     PAR(i,j,k) = 0.
801     do ilam = 1,tlam
802     PAR(i,j,k) = PAR(i,j,k) + PARw_k(ilam,k)
803     enddo
804     enddo
805    
806     #else /* DAR_RADTRANS */
807     c ------------ FULL RADIATIVE TRANSFER CODE ----------------------------
808     do ilam = 1,tlam
809     Edwsf(ilam) = oasim_ed(i,j,ilam,bi,bj)
810     Eswsf(ilam) = oasim_es(i,j,ilam,bi,bj)
811     enddo
812    
813     #ifdef DAR_RADTRANS_USE_MODEL_CALENDAR
814     C simplified solar zenith angle for 360-day year and daily averaged light
815     C cos(solz) is average over daylight period
816     call darwin_solz360(newtime, YC(i,j,bi,bj),
817     O solz)
818    
819     #else /* not DAR_RADTRANS_USE_MODEL_CALENDAR */
820     C use calendar date for full solar zenith angle computation
821     C Use local noon zenith angle to avoid problems with zero cosine and
822     C non-zero light. One should really use a zenith angle compatible with
823     C the light fields, in particular averaged over the same time period.
824     isec = MOD(36.*3600. - 240.*XC(i,j,bi,bj), 86400.)
825     call radtrans_sfcsolz(rad,iyr,imon,iday,isec,
826     I XC(i,j,bi,bj),YC(i,j,bi,bj),
827     O solz)
828     #endif /* not DAR_RADTRANS_USE_MODEL_CALENDAR */
829    
830     c have Ed,Es below surface - no need for this adjustment on Ed Es for surface affects
831     c do ilam=1,tlam
832     c rod(ilam) = 0.0 _d 0
833     c ros(ilam) = 0.0 _d 0
834     c enddo
835    
836     c compute 1/cos(zenith) for direct light below surface
837     call radtrans_sfcrmud(rad,solz,
838     O rmud)
839    
840     C compute absorption/scattering coefficients for radtrans
841     DO k=1,Nr
842     dz_k(k) = drF(k)*HFacC(i,j,k,bi,bj)
843     DO ilam = 1,tlam
844     c absorption by phyto
845     actot = 0.0
846     bctot = 0.0
847     bbctot = 0.0
848     DO np = 1,npmax
849     actot = actot + phychl_k(np,k)*aphy_chl(np,ilam)
850     bctot = bctot + phychl_k(np,k)*bphy_chl(np,ilam)
851     bbctot = bbctot + phychl_k(np,k)*bbphy_chl(np,ilam)
852     ENDDO
853     c particulate
854     apart_k(k,ilam) = part_k(k)*apart_P(ilam)
855     bpart_k(k,ilam) = part_k(k)*bpart_P(ilam)
856     bbpart_k(k,ilam) = part_k(k)*bbpart_P(ilam)
857     c add water and CDOM
858     a_k(k,ilam) = aw(ilam)+acdom_k(k,ilam)+actot+apart_k(k,ilam)
859     bt_k(k,ilam) = bw(ilam) + bctot + bpart_k(k,ilam)
860     bb_k(k,ilam) = darwin_bbw*bw(ilam)+bbctot+bbpart_k(k,ilam)
861     bb_k(k,ilam) = MAX(darwin_bbmin, bb_k(k,ilam))
862     ENDDO
863     ENDDO
864    
865     #ifdef DAR_RADTRANS_ITERATIVE
866     call darwin_radtrans_iter(
867     I dz_k,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
868     I darwin_radtrans_kmax,darwin_radtrans_niter,
869     O Edz,Esz,Euz,Eutop,
870     O tirrq,tirrwq,
871     I myThid)
872     #else
873     c dzlocal ?????
874     call darwin_radtrans(
875     I drF,rmud,Edwsf,Eswsf,a_k,bt_k,bb_k,
876     O Edz,Esz,Euz,Eutop,
877     O tirrq,tirrwq,
878     I myThid)
879     #endif
880     c
881     c uses chl from prev timestep (as wavebands does)
882     c keep like this in case need to consider upwelling irradiance as affecting the grid box above
883     c will pass to plankton: PARw only, but will be for this timestep for RT and prev timestep for WAVBANDS
884     c
885     c now copy
886     DO k=1,Nr
887     PAR(i,j,k) = tirrq(k)
888     DO ilam = 1,tlam
889     PARw_k(ilam,k) = tirrwq(ilam,k)
890     ENDDO
891     ENDDO
892     #endif /* DAR_RADTRANS */
893    
894     c oj: ???
895     c so PARw and PARwup from WAVEBANDS_1D are from previous timestep (attenuation done in plankton)
896     c but PARw and PARwup from WAVEBANDS_3D and RADTRANS are for the current timestep
897    
898     #endif /* WAVEBANDS */
899    
900     C ============================ k loop ==================================
901     c for each layer ...
902     do k= 1, NR
903     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
904    
905     c make sure we only deal with positive definite numbers
906     c brute force...
907     po4l = max(Ptr(i,j,k,bi,bj,iPO4 ),0. _d 0)
908     no3l = max(Ptr(i,j,k,bi,bj,iNO3 ),0. _d 0)
909     fetl = max(Ptr(i,j,k,bi,bj,iFeT ),0. _d 0)
910     sil = max(Ptr(i,j,k,bi,bj,iSi ),0. _d 0)
911     dopl = max(Ptr(i,j,k,bi,bj,iDOP ),0. _d 0)
912     donl = max(Ptr(i,j,k,bi,bj,iDON ),0. _d 0)
913     dofel = max(Ptr(i,j,k,bi,bj,iDOFe ),0. _d 0)
914     DO nz = 1,nzmax
915     ZooP(nz) = max(Ptr(i,j,k,bi,bj,iZooP (nz)),0. _d 0)
916     ZooN(nz) = max(Ptr(i,j,k,bi,bj,iZooN (nz)),0. _d 0)
917     ZooFe(nz) = max(Ptr(i,j,k,bi,bj,iZooFe(nz)),0. _d 0)
918     ZooSi(nz) = max(Ptr(i,j,k,bi,bj,iZooSi(nz)),0. _d 0)
919     ENDDO
920     popl = max(Ptr(i,j,k,bi,bj,iPOP ),0. _d 0)
921     ponl = max(Ptr(i,j,k,bi,bj,iPON ),0. _d 0)
922     pofel = max(Ptr(i,j,k,bi,bj,iPOFe ),0. _d 0)
923     psil = max(Ptr(i,j,k,bi,bj,iPOSi ),0. _d 0)
924     NH4l = max(Ptr(i,j,k,bi,bj,iNH4 ),0. _d 0)
925     NO2l = max(Ptr(i,j,k,bi,bj,iNO2 ),0. _d 0)
926     #ifdef ALLOW_CARBON
927     dicl = max(Ptr(i,j,k,bi,bj,iDIC ),0. _d 0)
928     docl = max(Ptr(i,j,k,bi,bj,iDOC ),0. _d 0)
929     pocl = max(Ptr(i,j,k,bi,bj,iPOC ),0. _d 0)
930     picl = max(Ptr(i,j,k,bi,bj,iPIC ),0. _d 0)
931     alkl = max(Ptr(i,j,k,bi,bj,iALK ),0. _d 0)
932     o2l = max(Ptr(i,j,k,bi,bj,iO2 ),0. _d 0)
933     cal = max(Ptr(i,j,k,bi,bj,iCa ),0. _d 0)
934    
935     DO nz = 1,nzmax
936     ZooCl(nz) = max(Ptr(i,j,k,bi,bj,iZooC (nz)),0. _d 0)
937     ENDDO
938     #endif
939    
940     totphyC = 0. _d 0
941     DO np=1,npmax
942     totphyC = totphyC + R_PC(np)*Ptr(i,j,k,bi,bj,iPhy+np-1)
943     ENDDO
944    
945     DO np = 1,npmax
946     Phy(np) = Phy_k(np,k)
947     #ifdef GEIDER
948     phychl(np) = phychl_k(np,k)
949     #endif
950     ENDDO
951    
952     #ifdef DAR_DIAG_DIVER
953     Diver1(i,j,k)=0. _d 0
954     Diver2(i,j,k)=0. _d 0
955     Diver3(i,j,k)=0. _d 0
956     Diver4(i,j,k)=0. _d 0
957     totphy=0. _d 0
958     do np=1,npmax
959     totphy=totphy + Phy(np)
960     tmpphy(np)=Phy(np)
961     enddo
962     if (totphy.gt.diver_thresh0) then
963     do np=1,npmax
964     c simple threshhold
965     if (Phy(np).gt.diver_thresh1) then
966     Diver1(i,j,k)=Diver1(i,j,k)+1. _d 0
967     endif
968     c proportion of total biomass
969     if (Phy(np)/totphy.gt.diver_thresh2) then
970     Diver2(i,j,k)=Diver2(i,j,k)+1. _d 0
971     endif
972     enddo
973     c majority of biomass by finding rank order
974     biotot=0. _d 0
975     do np2=1,npmax
976     phymax=0. _d 0
977     do np=1,npmax
978     if (tmpphy(np).gt.phymax) then
979     phymax=tmpphy(np)
980     npsave=np
981     endif
982     enddo
983     if (biotot.lt.totphy*diver_thresh3) then
984     Diver3(i,j,k)=Diver3(i,j,k)+1. _d 0
985     endif
986     biotot=biotot+tmpphy(npsave)
987     tmpphy(npsave)=0. _d 0
988     if (np2.eq.1) then
989     maxphy=phymax
990     endif
991     enddo
992     c ratio of maximum species
993     do np=1,npmax
994     if (Phy(np).gt.diver_thresh4*maxphy) then
995     Diver4(i,j,k)=Diver4(i,j,k)+1. _d 0
996     endif
997     enddo
998     endif
999     #endif
1000    
1001     c..........................................................
1002     c find local light
1003     c..........................................................
1004    
1005     PARl = PAR(i,j,k)
1006     c..........................................................
1007    
1008     c for explicit sinking of particulate matter and phytoplankton
1009     if (k.eq.1) then
1010     popupl =0. _d 0
1011     ponupl =0. _d 0
1012     pofeupl = 0. _d 0
1013     psiupl = 0. _d 0
1014     do np=1,npmax
1015     Phyup(np)=0. _d 0
1016     #ifdef DYNAMIC_CHL
1017     chlup(np)=0. _d 0
1018     #endif
1019     enddo
1020     #ifdef ALLOW_CARBON
1021     pocupl = 0. _d 0
1022     picupl = 0. _d 0
1023     #endif
1024     endif
1025    
1026     #ifdef ALLOW_LONGSTEP
1027     Tlocal = LS_theta(i,j,k,bi,bj)
1028     Slocal = LS_salt(i,j,k,bi,bj)
1029     #else
1030     Tlocal = theta(i,j,k,bi,bj)
1031     Slocal = salt(i,j,k,bi,bj)
1032     #endif
1033    
1034     freefu = max(freefe(i,j,k),0. _d 0)
1035     if (k.eq.1) then
1036     inputFel = inputFe(i,j,bi,bj)
1037     else
1038     inputFel = 0. _d 0
1039     endif
1040    
1041     dzlocal = drF(k)*HFacC(i,j,k,bi,bj)
1042     c set bottom=1.0 if the layer below is not ocean
1043     ktmp=min(nR,k+1)
1044     if(hFacC(i,j,ktmp,bi,bj).eq.0. _d 0.or.k.eq.Nr) then
1045     bottom = 1.0 _d 0
1046     else
1047     bottom = 0.0 _d 0
1048     endif
1049    
1050     c set tendencies to 0
1051     do np=1,npmax
1052     dphy(np)=0. _d 0
1053     enddo
1054     do nz=1,nzmax
1055     dzoop(nz)=0. _d 0
1056     dzoon(nz)=0. _d 0
1057     dzoofe(nz)=0. _d 0
1058     dzoosi(nz)=0. _d 0
1059     enddo
1060     dPO4l=0. _d 0
1061     dNO3l=0. _d 0
1062     dFeTl=0. _d 0
1063     dSil=0. _d 0
1064     dDOPl=0. _d 0
1065     dDONl=0. _d 0
1066     dDOFel=0. _d 0
1067     dPOPl=0. _d 0
1068     dPONl=0. _d 0
1069     dPOFel=0. _d 0
1070     dPSil=0. _d 0
1071     dNH4l=0. _d 0
1072     dNO2l=0. _d 0
1073     #ifdef DYNAMIC_CHL
1074     do np=1,npmax
1075     dphychl(np)=0. _d 0
1076     enddo
1077     #endif
1078     #ifdef ALLOW_CARBON
1079     ddicl=0. _d 0
1080     ddocl=0. _d 0
1081     dpocl=0. _d 0
1082     dpicl=0. _d 0
1083     dalkl=0. _d 0
1084     do2l=0. _d 0
1085     do nz=1,nzmax
1086     dzoocl(nz)=0. _d 0
1087     enddo
1088     #endif
1089     #ifdef ADKINS_SURF_FLUX
1090     dcal = 0. _d 0
1091     #endif
1092     c set other arguments to zero
1093     PP=0. _d 0
1094     Nfix=0. _d 0
1095     denit=0. _d 0
1096     do np=1,npmax
1097     Rstarl(np)=0. _d 0
1098     RNstarl(np)=0. _d 0
1099     #ifdef DAR_DIAG_GROW
1100     Growl(np)=0. _d 0
1101     Growsql(np)=0. _d 0
1102     #endif
1103     #ifdef ALLOW_DIAZ
1104     #ifdef DAR_DIAG_NFIXP
1105     NfixPl(np)=0. _d 0
1106     #endif
1107     #endif
1108     enddo
1109    
1110     c if (i.eq.20.and.j.eq.20.and.k.eq.1) debug=8
1111     c if (i.eq.10.and.j.eq.10.and.k.eq.1) debug=100
1112     c if (i.eq.1.and.j.eq.10.and.k.eq.1) debug=10
1113     c if (i.eq.1.and.j.eq.1.and.k.eq.10) debug=14
1114    
1115     if (debug.eq.7) print*,'PO4, DOP, POP, ZooP',
1116     & PO4l, DOPl, POPl, zooP
1117     if (debug.eq.7) print*,'NO3, NO2, NH4, DON, PON, ZooN',
1118     & NO3l,NO2l,NH4l, DONl, PONl, ZooN
1119     if (debug.eq.7) print*,'FeT, DOFe, POFe, Zoofe',
1120     & FeTl, DOFel, POFel, zooFe
1121     if (debug.eq.7) print*,'Si, Psi, zooSi',
1122     & Sil, PSil, zooSi
1123     if (debug.eq.7) print*,'Total Phy', sumpy, PARl, lite
1124     if (debug.eq.7) print*,'Phy', Phy
1125    
1126     if (debug.eq.8) print*,'k, PARl, inputFel, dzlocal',
1127     & PARl, inputFel, dzlocal
1128    
1129     c if (NO3l.eq.0. _d 0.or.NO2l.eq.0. _d 0
1130     c & .or.NH4l.eq.0. _d 0) then
1131     c print*,'QQ N zeros',i,j,k,NO3l,NO2l,NH4l
1132     c endif
1133    
1134     c compute Ksp as a function of temperature and pressure
1135    
1136     bdepth = 0.0d0
1137     cdepth = 0.0d0
1138     pressc = 1.0d0
1139    
1140     do l = 1,k
1141     cdepth = bdepth + 0.5d0*drF(l)
1142     bdepth = bdepth + drF(l)
1143     pressc = 1.0d0 + 0.1d0*cdepth
1144     end do
1145    
1146     if (maskC(i,j,k,bi,bj).NE.0. _d 0) then
1147     t = Tlocal
1148     s = max(4. _d 0, Slocal)
1149    
1150     tk = 273.15 + t
1151     tk100 = tk/100.0
1152     tk1002=tk100*tk100
1153     invtk=1.0/tk
1154     dlogtk=log(tk)
1155     is=19.924*s/(1000.-1.005*s)
1156     is2=is*is
1157     sqrtis=sqrt(is)
1158     s2=s*s
1159     sqrts=sqrt(s)
1160     s15=s**1.5
1161     scl=s/1.80655
1162    
1163     c f = k0(1-pH2O)*correction term for non-ideality
1164     c Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
1165     ff(i,j,bi,bj) = exp(-162.8301 + 218.2968/tk100 +
1166     & 90.9241*log(tk100)-1.47696*tk1002 +
1167     & s*(0.025695-0.025225*tk100 +
1168     & 0.0049867*tk1002))
1169    
1170     c K0 from Weiss 1974
1171     ak0(i,j,bi,bj) = exp(93.4517/tk100-60.2409 +
1172     & 23.3585*log(tk100) +
1173     & s*(0.023517-0.023656*tk100 +
1174     & 0.0047036*tk1002))
1175    
1176     c k1 = [H][HCO3]/[H2CO3]
1177     c k2 = [H][CO3]/[HCO3]
1178     c Millero p.664 (1995) using Mehrbach et al. data on seawater scale
1179     ak1(i,j,bi,bj)=10**(-1*(3670.7*invtk -
1180     & 62.008+9.7944*dlogtk -
1181     & 0.0118*s+0.000116*s2))
1182    
1183     ak2(i,j,bi,bj)=10**(-1*(1394.7*invtk+4.777 -
1184     & 0.0184*s+0.000118*s2))
1185    
1186     c NOW PRESSURE DEPENDENCE:
1187     c Following Takahashi (1981) GEOSECS report - quoting Culberson and
1188     c Pytkowicz (1968)
1189     c pressc = pressure in bars
1190     ak1(i,j,bi,bj) = ak1(i,j,bi,bj) *
1191     & exp((24.2-0.085*t)*(pressc-1.0)/(83.143*tk))
1192    
1193     c FIRST GO FOR K2: According to GEOSECS (1982) report
1194     c ak2(i,j,bi,bj) = ak2(i,j,bi,bj) *
1195     c & exp((26.4-0.040*t)*(pressc-1.0)/(83.143*tk))
1196    
1197     c SECOND GO FOR K2: corrected coeff according to CO2sys documentation
1198     c E. Lewis and D. Wallace (1998) ORNL/CDIAC-105
1199     ak2(i,j,bi,bj) = ak2(i,j,bi,bj) *
1200     & exp((16.4-0.040*t)*(pressc-1.0)/(83.143*tk))
1201    
1202     c kb = [H][BO2]/[HBO2]
1203     c Millero p.669 (1995) using data from dickson (1990)
1204     akb(i,j,bi,bj)=exp((-8966.90-2890.53*sqrts-77.942*s +
1205     & 1.728*s15-0.0996*s2)*invtk +
1206     & (148.0248+137.1942*sqrts+1.62142*s) +
1207     & (-24.4344-25.085*sqrts-0.2474*s) *
1208     & dlogtk+0.053105*sqrts*tk)
1209    
1210     c Mick and Karsten - Dec 04
1211     c ADDING pressure dependence based on Millero (1995), p675
1212     c with additional info from CO2sys documentation (E. Lewis and
1213     c D. Wallace, 1998 - see endnotes for commentary on Millero, 95)
1214     bigR = 83.145
1215     dv = -29.48+0.1622*t+2.608d-3*t*t
1216     dk = -2.84d-3
1217     pfactor = - (dv/(bigR*tk))*pressc
1218     & + (0.5*dk/(bigR*tk))*pressc*pressc
1219    
1220     akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)
1221    
1222     c k1p = [H][H2PO4]/[H3PO4]
1223     c DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
1224     ak1p(i,j,bi,bj) = exp(-4576.752*invtk+115.525 -
1225     & 18.453*dlogtk +
1226     & (-106.736*invtk+0.69171)*sqrts +
1227     & (-0.65643*invtk-0.01844)*s)
1228    
1229     c k2p = [H][HPO4]/[H2PO4]
1230     c DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
1231     ak2p(i,j,bi,bj) = exp(-8814.715*invtk+172.0883 -
1232     & 27.927*dlogtk +
1233     & (-160.340*invtk+1.3566)*sqrts +
1234     & (0.37335*invtk-0.05778)*s)
1235    
1236     c k3p = [H][PO4]/[HPO4]
1237     c DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
1238     ak3p(i,j,bi,bj) = exp(-3070.75*invtk-18.141 +
1239     & (17.27039*invtk+2.81197) *
1240     & sqrts+(-44.99486*invtk-0.09984)*s)
1241    
1242     c ksi = [H][SiO(OH)3]/[Si(OH)4]
1243     c Millero p.671 (1995) using data from Yao and Millero (1995)
1244     aksi(i,j,bi,bj) = exp(-8904.2*invtk+117.385 -
1245     & 19.334*dlogtk +
1246     & (-458.79*invtk+3.5913)*sqrtis +
1247     & (188.74*invtk-1.5998)*is +
1248     & (-12.1652*invtk+0.07871)*is2 +
1249     & log(1.0-0.001005*s))
1250    
1251     c kw = [H][OH]
1252     c Millero p.670 (1995) using composite data
1253     akw(i,j,bi,bj) = exp(-13847.26*invtk+148.9652 -
1254     & 23.6521*dlogtk +
1255     & (118.67*invtk-5.977+1.0495*dlogtk) *
1256     & sqrts-0.01615*s)
1257    
1258     c ks = [H][SO4]/[HSO4]
1259     c dickson (1990, J. chem. Thermodynamics 22, 113)
1260     aks(i,j,bi,bj)=exp(-4276.1*invtk+141.328 -
1261     & 23.093*dlogtk +
1262     & (-13856*invtk+324.57-47.986*dlogtk)*sqrtis +
1263     & (35474*invtk-771.54+114.723*dlogtk)*is -
1264     & 2698*invtk*is**1.5+1776*invtk*is2 +
1265     & log(1.0-0.001005*s))
1266    
1267     c kf = [H][F]/[HF]
1268     c dickson and Riley (1979) -- change pH scale to total
1269     akf(i,j,bi,bj)=exp(1590.2*invtk-12.641+1.525*sqrtis +
1270     & log(1.0-0.001005*s) +
1271     & log(1.0+(0.1400/96.062)*(scl)/aks(i,j,bi,bj)))
1272    
1273     c Calculate concentrations for borate, sulfate, and fluoride
1274     c Uppstrom (1974)
1275     bt(i,j,bi,bj) = 0.000232*scl/10.811
1276    
1277     c Morris & Riley (1966)
1278     st(i,j,bi,bj) = 0.14*scl/96.062
1279    
1280     c Riley (1965)
1281     ft(i,j,bi,bj) = 0.000067*scl/18.9984
1282    
1283     c solubility product for calcite
1284     C Following Mucci (1983) - from Zeebe/Wolf-Gladrow equic.m
1285     tmpa1 = -171.9065-(0.077993*tk)+(2839.319/tk) +
1286     & (71.595*log10(tk))
1287    
1288     tmpa2 = +(-0.77712+(0.0028426*tk)+(178.34/tk))*sqrts
1289     tmpa3 = -(0.07711*s)+(0.0041249*s15)
1290     logKspc = tmpa1+tmpa2+tmpa3
1291     Ksp_T_Calc = 10.0**logKspc
1292    
1293     c alternative pressure dependence from Ingle (1975)
1294     zdum = (pressc*10.0d0-10.0d0)/10.0d0
1295     xvalue = ((48.8d0-0.53d0*t)*zdum +
1296     & (-0.00588d0+0.0001845d0*t)*zdum*zdum) /
1297     & (188.93d0*(t+273.15d0))
1298    
1299     KspTP(i,j,bi,bj) = Ksp_T_Calc*10**(xvalue)
1300    
1301     else
1302     ff(i,j,bi,bj)=0.d0
1303     ak0(i,j,bi,bj)= 0.d0
1304     ak1(i,j,bi,bj)= 0.d0
1305     ak2(i,j,bi,bj)= 0.d0
1306     akb(i,j,bi,bj)= 0.d0
1307     ak1p(i,j,bi,bj) = 0.d0
1308     ak2p(i,j,bi,bj) = 0.d0
1309     ak3p(i,j,bi,bj) = 0.d0
1310     aksi(i,j,bi,bj) = 0.d0
1311     akw(i,j,bi,bj) = 0.d0
1312     aks(i,j,bi,bj)= 0.d0
1313     akf(i,j,bi,bj)= 0.d0
1314     bt(i,j,bi,bj) = 0.d0
1315     st(i,j,bi,bj) = 0.d0
1316     ft(i,j,bi,bj) = 0.d0
1317     KspTP(i,j,bi,bj) = 0.d0
1318     endif
1319    
1320     c compute CO3 for DARWIN_PLANKTON and sediment fluxes
1321     if (maskC(i,j,k,bi,bj).NE.0. _d 0) then
1322    
1323     pCO2SolverTemp = max(-4. _d 0, min(39. _d 0, Tlocal))
1324     pCO2SolverSal = max(4. _d 0, min(50. _d 0, Slocal))
1325    
1326     c check bounds for PCO2 solver
1327     pCO2SolverDic = max(100. _d 0, min(4000. _d 0,dicl))
1328     pCO2SolverAlk = max(100. _d 0, min(4000. _d 0,alkl))
1329     pCO2SolverPo4 = max(1. _d -10, min(10. _d 0,po4l))
1330     pCO2SolverSi = max(1. _d -8, min(500. _d 0,sil))
1331    
1332     c convert to mol m^-3
1333     pCO2SolverDic = pCO2SolverDic*1.0 _d -3
1334     pCO2SolverAlk = pCO2SolverAlk*1.0 _d -3
1335     pCO2SolverPo4 = pCO2SolverPo4*1.0 _d -3
1336     pCO2SolverSi = pCO2SolverSi*1.0 _d -3
1337    
1338     CALL CALC_PCO2_APPROX(
1339     I pCO2SolverTemp,pCO2SolverSal,
1340     I pCO2SolverDic,pCO2SolverPo4,
1341     I pCO2SolverSi,pCO2SolverAlk,
1342     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
1343     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
1344     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
1345     I aksi(i,j,bi,bj),akf(i,j,bi,bj),
1346     I ak0(i,j,bi,bj),fugf(i,j,bi,bj),
1347     I ff(i,j,bi,bj),
1348     I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
1349     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),CO3(i,j,bi,bj),
1350     I myThid)
1351    
1352     else
1353    
1354     pH(i,j,bi,bj) = 0. _d 0
1355     pCO2(i,j,bi,bj) = 0. _d 0
1356     CO3(i,j,bi,bj) = 0. _d 0
1357    
1358     endif
1359    
1360     #ifdef ADKINS_SURF_FLUX
1361     c include surface fluxes from forcing files
1362     if (k.eq.1) then
1363     dcal = dcal + (caSurf_flx(i,j,bi,bj) * 1. _d 3)
1364     dcal = dcal - budgetConsumpDIC_PIC(i,j,k,bi,bj) +
1365     & disscPIC(i,j,k,bi,bj)
1366     endif
1367     #endif /* ADKINS_SURF_FLUX */
1368    
1369     Ptr(i,j,k,bi,bj,iCA) = Ptr(i,j,k,bi,bj,iCA) +
1370     & dtplankton*dcal
1371    
1372     calcium(i,j,bi,bj) = Ptr(i,j,k,bi,bj,iCa)
1373    
1374     c ANNA pass extra variables if WAVEBANDS
1375     CALL DARWIN_PLANKTON(
1376     U Phy,
1377     I zooP, zooN, zooFe, zooSi,
1378     O PP, Chl, Nfix, denit,
1379     I PO4l, NO3l, FeTl, Sil,
1380     I NO2l, NH4l,
1381     I DOPl, DONl, DOFel,
1382     I POPl, PONl, POFel, PSil,
1383     I phyup, popupl, ponupl,
1384     I pofeupl, psiupl,
1385     I PARl,
1386     I Tlocal, Slocal,
1387     I freefu, inputFel,
1388     I bottom, dzlocal,
1389     O Rstarl, RNstarl,
1390     #ifdef DAR_DIAG_GROW
1391     O Growl, Growsql,
1392     #endif
1393     #ifdef ALLOW_DIAZ
1394     #ifdef DAR_DIAG_NFIXP
1395     O NfixPl,
1396     #endif
1397     #endif
1398     O dphy, dzooP, dzooN, dzooFe,
1399     O dzooSi,
1400     O dPO4l, dNO3l, dFeTl, dSil,
1401     O dNH4l, dNO2l,
1402     O dDOPl, dDONl, dDOFel,
1403     O dPOPl, dPONl, dPOFel, dPSil,
1404     #ifdef ALLOW_CARBON
1405     I dicl, docl, pocl, picl,
1406     I alkl, o2l, zoocl,
1407     I pocupl, picupl, KspTP(i,j,bi,bj),
1408     I CO3(i,j,bi,bj),calcium(i,j,bi,bj),
1409     O ddicl, ddocl, dpocl, dpicl,
1410     O dalkl, do2l, dzoocl,omegaC(i,j,k,bi,bj),
1411     O disscPIC(i,j,k,bi,bj),
1412     #endif
1413     #ifdef CO2_FLUX_BUDGET
1414     O budgetConsumpDIC(i,j,k,bi,bj),
1415     O budgetConsumpDIC_PIC(i,j,k,bi,bj),
1416     O budgetDOCRemin(i,j,k,bi,bj),
1417     O budgetPReminC(i,j,k,bi,bj),
1418     #endif /* CO2_FLUX_BUDGET */
1419     #ifdef GEIDERCO3
1420     O phychl,
1421     #ifdef DYNAMIC_CHL
1422     I dphychl,
1423     I chlup,
1424     #endif
1425     #ifdef WAVEBANDS
1426     I PARw_k(1,k),
1427     #endif
1428     #endif
1429     #ifdef ALLOW_PAR_DAY
1430     I PARday(i,j,k,bi,bj,PARiprev),
1431     #endif
1432     #ifdef DAR_DIAG_CHL
1433     O ChlGeiderlocal, ChlDoneylocal,
1434     O ChlCloernlocal,
1435     #endif
1436     I debug,
1437     I runtim,
1438     I MyThid)
1439    
1440     #ifdef IRON_SED_SOURCE
1441     c only above minimum depth (continental shelf)
1442     if (rF(k).gt.-depthfesed) then
1443     c only if bottom layer
1444     if (bottom.eq.1.0 _d 0) then
1445     #ifdef IRON_SED_SOURCE_VARIABLE
1446     c calculate sink of POP into bottom layer
1447     tmp=(wp_sink*POPupl)/(dzlocal)
1448     c convert to dPOCl
1449     dFetl=dFetl+fesedflux_pcm*(tmp*106. _d 0)
1450     #else
1451     dFetl=dFetl+fesedflux/
1452     & (drF(k)*hFacC(i,j,k,bi,bj))
1453     #endif
1454     endif
1455     endif
1456     #endif
1457     ponupl = PONl
1458     pofeupl = POFel
1459     psiupl = PSil
1460     do np=1,npmax
1461     Phyup(np) = Phy(np)
1462     #ifdef DYNAMIC_CHL
1463     chlup(np) = phychl(np)
1464     #endif
1465     enddo
1466    
1467     c
1468     #ifdef ALLOW_CARBON
1469     pocupl = POCl
1470     picupl = PICl
1471     c include surface forcing
1472     if (k.eq.1) then
1473     ddicl = ddicl + flxCO2(i,j)
1474     dalkl = dalkl + flxALK(i,j)
1475     do2l = do2l + flxO2(i,j)
1476     #ifdef ADKINS_SURF_FLUX
1477     c include surface fluxes from forcing files
1478     ddicl = ddicl + (dicSurf_flx(i,j,bi,bj) * 1. _d 3)
1479     dalkl = dalkl + (alkSurf_flx(i,j,bi,bj) * 1. _d 3)
1480     #endif /* ADKINS_SURF_FLUX */
1481     endif
1482     #ifdef ALLOW_SED_DISS_FLUX
1483     c compute sediment dissolution fluxes
1484     if (bottom.eq.1.0 _d 0) then
1485    
1486     CO3Sed(i,j,bi,bj) = (KspTP(i,j,bi,bj) /
1487     & calcium(i,j,bi,bj))
1488    
1489     BBLDiffusionCoeffLoc = 4.01 _d -10
1490    
1491     if (darwin_BBLFile .NE. ' ') then
1492     BBLThicknessLoc = BBLThickness(i,j,bi,bj)
1493     else
1494     BBLThicknessLoc = 1.0 _d -3
1495     endif
1496    
1497     CO3Sw(i,j,bi,bj) = CO3(i,j,bi,bj)
1498    
1499     c print*,'CO3Sw: ',CO3Sw(i,j,bi,bj)
1500    
1501     c convert from mol kg^-1 to mol m^-3
1502     CO3Sed(i,j,bi,bj) = CO3Sed(i,j,bi,bj) * rhoConst
1503     CO3Sw(i,j,bi,bj) = CO3Sw(i,j,bi,bj) * rhoConst
1504    
1505     c compute flux in mol m^-3 s^-1
1506     RFlux(i,j,bi,bj) = (BBLDiffusionCoeffLoc /
1507     & BBLThicknessLoc) * (CO3Sed(i,j,bi,bj)-CO3Sw(i,j,bi,bj))
1508    
1509     c prevent negative fluxes (for now)
1510     if(RFlux(i,j,bi,bj).LT.0) then
1511     RFlux(i,j,bi,bj) = 0;
1512     endif
1513    
1514     DICSedFlux(i,j,bi,bj) = RFlux(i,j,bi,bj)
1515     ALKSedFlux(i,j,bi,bj) = 2*RFlux(i,j,bi,bj)
1516    
1517     c convert sediment fluxes to mmol m^-3 s^-1 and then add to local values
1518     ddicl = ddicl+(RFlux(i,j,bi,bj)*1.0 _d 3 /
1519     & drF(k)*HFacC(i,j,k,bi,bj))
1520     dalkl = dalkl+(2*RFlux(i,j,bi,bj)*1.0 _d 3 /
1521     & drF(k)*HFacC(i,j,k,bi,bj))
1522    
1523     endif
1524     #endif /* ALLOW_SED_DISS_FLUX */
1525     #endif /* ALLOW_CARBON */
1526    
1527     #ifdef CONS_SUPP
1528     c only works for two layer model
1529     if (k.eq.2) then
1530     dpo4l=0. _d 0
1531     dno3l=0. _d 0
1532     dfetl=0. _d 0
1533     dsil=0. _d 0
1534     endif
1535     #endif
1536     #ifdef RELAX_NUTS
1537     #ifdef DENIT_RELAX
1538     if (rF(k).lt.-depthdenit) then
1539     if (darwin_relaxscale.gt.0. _d 0) then
1540     IF ( darwin_NO3_RelaxFile .NE. ' ' ) THEN
1541     c Fanny's formulation
1542     tmp=(Ptr(i,j,k,bi,bj,iNO3 )-no3_obs(i,j,k,bi,bj))
1543     if (tmp.gt.0. _d 0) then
1544     dno3l=dno3l-(tmp/
1545     & darwin_relaxscale)
1546     denit=tmp/
1547     & darwin_relaxscale
1548     else
1549     denit=0. _d 0
1550     endif
1551     c --- end fanny's formulation
1552     ENDIF
1553     c steph's alternative
1554     c tmp=(Ptr(i,j,k,bi,bj,iNO3 )-
1555     c & 16. _d 0 * Ptr(i,j,k,bi,bj,iPO4 ))
1556     c if (tmp.gt.0. _d 0) then
1557     c dno3l=dno3l-(tmp/
1558     c & darwin_relaxscale)
1559     c denit=tmp/
1560     c & darwin_relaxscale
1561     c else
1562     c denit=0. _d 0
1563     c endif
1564     c ---- end steph's alternative
1565     endif
1566     endif
1567     #else
1568     if (darwin_relaxscale.gt.0. _d 0) then
1569     IF ( darwin_PO4_RelaxFile .NE. ' ' ) THEN
1570     tmp=(Ptr(i,j,k,bi,bj,iPO4 )-po4_obs(i,j,k,bi,bj))
1571     if (tmp.lt.0. _d 0) then
1572     dpo4l=dpo4l-(tmp/
1573     & darwin_relaxscale)
1574     endif
1575     ENDIF
1576     IF ( darwin_NO3_RelaxFile .NE. ' ' ) THEN
1577     tmp=(Ptr(i,j,k,bi,bj,iNO3 )-no3_obs(i,j,k,bi,bj))
1578     if (tmp.lt.0. _d 0) then
1579     dno3l=dno3l-(tmp/
1580     & darwin_relaxscale)
1581     endif
1582     ENDIF
1583     IF ( darwin_Fet_RelaxFile .NE. ' ' ) THEN
1584     tmp=(Ptr(i,j,k,bi,bj,iFeT )-fet_obs(i,j,k,bi,bj))
1585     if (tmp.lt.0. _d 0) then
1586     dfetl=dfetl-(tmp/
1587     & darwin_relaxscale)
1588     endif
1589     ENDIF
1590     IF ( darwin_Si_RelaxFile .NE. ' ' ) THEN
1591     tmp=( Ptr(i,j,k,bi,bj,iSi )-si_obs(i,j,k,bi,bj))
1592     if (tmp.lt.0. _d 0) then
1593     dsil=dsil-(tmp/
1594     & darwin_relaxscale)
1595     endif
1596     ENDIF
1597     endif
1598     #endif
1599     #endif
1600     #ifdef FLUX_NUTS
1601     dpo4l=dpo4l+po4_flx(i,j,k,bi,bj)
1602     dno3l=dno3l+no3_flx(i,j,k,bi,bj)
1603     dfetl=dfetl+fet_flx(i,j,k,bi,bj)
1604     dsil=dsil+si_flx(i,j,k,bi,bj)
1605     #endif
1606    
1607     #ifdef ALLOW_OBCS
1608     IF (useOBCS) THEN
1609     dpo4l = dpo4l *maskInC(i,j,bi,bj)
1610     dno3l = dno3l *maskInC(i,j,bi,bj)
1611     dfetl = dfetl *maskInC(i,j,bi,bj)
1612     dsil = dsil *maskInC(i,j,bi,bj)
1613     ddopl = ddopl *maskInC(i,j,bi,bj)
1614     ddonl = ddonl *maskInC(i,j,bi,bj)
1615     ddofel = ddofel*maskInC(i,j,bi,bj)
1616     dpopl = dpopl *maskInC(i,j,bi,bj)
1617     dponl = dponl *maskInC(i,j,bi,bj)
1618     dpofel = dpofel*maskInC(i,j,bi,bj)
1619     dpsil = dpsil *maskInC(i,j,bi,bj)
1620     dnh4l = dnh4l *maskInC(i,j,bi,bj)
1621     dno2l = dno2l *maskInC(i,j,bi,bj)
1622     DO nz = 1,nzmax
1623     dzoop (nz) = dzoop (nz)*maskInC(i,j,bi,bj)
1624     dzoon (nz) = dzoon (nz)*maskInC(i,j,bi,bj)
1625     dzoofe(nz) = dzoofe(nz)*maskInC(i,j,bi,bj)
1626     dzoosi(nz) = dzoosi(nz)*maskInC(i,j,bi,bj)
1627     ENDDO
1628     DO np = 1,npmax
1629     dPhy(np) = dPhy(np)*maskInC(i,j,bi,bj)
1630     #ifdef GEIDER
1631     #ifdef DYNAMIC_CHL
1632     dphychl(np) = dphychl(np)*maskInC(i,j,bi,bj)
1633     #endif
1634     #endif
1635     ENDDO
1636     #ifdef ALLOW_CARBON
1637     ddicl = ddicl*maskInC(i,j,bi,bj)
1638     ddocl = ddocl*maskInC(i,j,bi,bj)
1639     dpocl = dpocl*maskInC(i,j,bi,bj)
1640     dpicl = dpicl*maskInC(i,j,bi,bj)
1641     dalkl = dalkl*maskInC(i,j,bi,bj)
1642     do2l = do2l *maskInC(i,j,bi,bj)
1643     dcal = dcal *maskInC(i,j,bi,bj)
1644     DO nz = 1,nzmax
1645     dzoocl(nz) = dzoocl(nz)*maskInC(i,j,bi,bj)
1646     ENDDO
1647     #endif
1648     ENDIF
1649     #endif
1650    
1651     c now update main tracer arrays
1652     dtplankton = PTRACERS_dTLev(k)/float(nsubtime)
1653     Ptr(i,j,k,bi,bj,iPO4 ) = Ptr(i,j,k,bi,bj,iPO4) +
1654     & dtplankton*dpo4l
1655     Ptr(i,j,k,bi,bj,iNO3 ) = Ptr(i,j,k,bi,bj,iNO3) +
1656     & dtplankton*dno3l
1657     Ptr(i,j,k,bi,bj,iFeT ) = Ptr(i,j,k,bi,bj,iFeT) +
1658     & dtplankton*dfetl
1659     Ptr(i,j,k,bi,bj,iSi ) = Ptr(i,j,k,bi,bj,iSi ) +
1660     & dtplankton*dsil
1661     Ptr(i,j,k,bi,bj,iDOP ) = Ptr(i,j,k,bi,bj,iDOP) +
1662     & dtplankton*ddopl
1663     Ptr(i,j,k,bi,bj,iDON ) = Ptr(i,j,k,bi,bj,iDON) +
1664     & dtplankton*ddonl
1665     Ptr(i,j,k,bi,bj,iDOFe) = Ptr(i,j,k,bi,bj,iDOFe) +
1666     & dtplankton*ddofel
1667     Ptr(i,j,k,bi,bj,iPOP ) = Ptr(i,j,k,bi,bj,iPOP ) +
1668     & dtplankton*dpopl
1669     Ptr(i,j,k,bi,bj,iPON ) = Ptr(i,j,k,bi,bj,iPON ) +
1670     & dtplankton*dponl
1671     Ptr(i,j,k,bi,bj,iPOFe) = Ptr(i,j,k,bi,bj,iPOFe) +
1672     & dtplankton*dpofel
1673     Ptr(i,j,k,bi,bj,iPOSi) = Ptr(i,j,k,bi,bj,iPOSi) +
1674     & dtplankton*dpsil
1675     Ptr(i,j,k,bi,bj,iNH4 ) = Ptr(i,j,k,bi,bj,iNH4 ) +
1676     & dtplankton*dnh4l
1677     Ptr(i,j,k,bi,bj,iNO2 ) = Ptr(i,j,k,bi,bj,iNO2 ) +
1678     & dtplankton*dno2l
1679     DO nz = 1,nzmax
1680     Ptr(i,j,k,bi,bj,iZooP (nz)) = Ptr(i,j,k,bi,bj,iZooP (nz)) +
1681     & dtplankton*dzoop (nz)
1682     Ptr(i,j,k,bi,bj,iZooN (nz)) = Ptr(i,j,k,bi,bj,iZooN (nz)) +
1683     & dtplankton*dzoon (nz)
1684     Ptr(i,j,k,bi,bj,iZooFe(nz)) = Ptr(i,j,k,bi,bj,iZooFe(nz)) +
1685     & dtplankton*dzoofe(nz)
1686     Ptr(i,j,k,bi,bj,iZooSi(nz)) = Ptr(i,j,k,bi,bj,iZooSi(nz)) +
1687     & dtplankton*dzoosi(nz)
1688     ENDDO
1689     DO np = 1,npmax
1690     Ptr(i,j,k,bi,bj,iPhy+np-1) = Ptr(i,j,k,bi,bj,iPhy+np-1) +
1691     & dtplankton*dPhy(np)
1692     #ifdef GEIDER
1693     #ifdef DYNAMIC_CHL
1694     if (np.eq.1) Chl=0. _d 0
1695     Ptr(i,j,k,bi,bj,iChl+np-1) = Ptr(i,j,k,bi,bj,iChl+np-1) +
1696     & dtplankton*dphychl(np)
1697     c chltmp=Ptr(i,j,k,bi,bj,iChl+np-1)
1698     c phytmp=Ptr(i,j,k,bi,bj,iPhy+np-1)
1699     c Ptr(i,j,k,bi,bj,iChl+np-1)=
1700     c & max(chltmp,phytmp*R_PC(np)*chl2cmin(np))
1701     c if (np.eq.1.and.i.eq.1.and.j.eq.1.and.k.eq.1)
1702     c & print*,chltmp,phytmp,phytmp*R_PC(np)*chl2cmin(np),
1703     c & phytmp*R_PC(np)*chl2cmax(np)
1704     c in darwin_plankton this is stored for previous timestep. Reset here.
1705     Chl=Chl+Ptr(i,j,k,bi,bj,iChl+np-1)
1706     #else
1707     Chl_phy(i,j,k,bi,bj,np)=phychl(np)
1708     #endif
1709     #endif
1710     ENDDO
1711     #ifdef ALLOW_CARBON
1712     Ptr(i,j,k,bi,bj,iDIC ) = Ptr(i,j,k,bi,bj,iDIC ) +
1713     & dtplankton*ddicl
1714     Ptr(i,j,k,bi,bj,iDOC ) = Ptr(i,j,k,bi,bj,iDOC ) +
1715     & dtplankton*ddocl
1716     Ptr(i,j,k,bi,bj,iPOC ) = Ptr(i,j,k,bi,bj,iPOC ) +
1717     & dtplankton*dpocl
1718     Ptr(i,j,k,bi,bj,iPIC ) = Ptr(i,j,k,bi,bj,iPIC ) +
1719     & dtplankton*dpicl
1720     Ptr(i,j,k,bi,bj,iALK ) = Ptr(i,j,k,bi,bj,iALK ) +
1721     & dtplankton*dalkl
1722     Ptr(i,j,k,bi,bj,iO2 ) = Ptr(i,j,k,bi,bj,iO2 ) +
1723     & dtplankton*do2l
1724     DO nz = 1,nzmax
1725     Ptr(i,j,k,bi,bj,iZooC (nz)) = Ptr(i,j,k,bi,bj,iZooC (nz)) +
1726     & dtplankton*dzoocl (nz)
1727     ENDDO
1728     #endif
1729     c
1730     #ifdef ALLOW_MUTANTS
1731     cQQQQTEST
1732     if (debug.eq.11) then
1733     if (k.lt.8) then
1734     do np=1,60
1735     if(mod(np,4).eq. 1. _d 0)then
1736     np2=np+1
1737     np4=np+3
1738    
1739     Coj: couldn't test this part after change Phynp -> Ptr(...,iPhy+np-1)
1740     Coj: used to be many copies of this:
1741     C if (dPhy(2).gt.dPhy(4).and.dPhy(4).gt.0. _d 0) then
1742     C print*,'QQQ dphy2 > dphy4',i,j,k,Phy2(i,j,k),
1743     C & Phy4(i,j,k), dPhy(2), dPhy(4)
1744     C endif
1745     C if (Phy2(i,j,k).gt.Phy4(i,j,k).and.
1746     C & Phy4(i,j,k).gt.0. _d 0) then
1747     C print*,'QQ phy02 > phy04',i,j,k,Phy2(i,j,k),
1748     C & Phy4(i,j,k), dPhy(2), dPhy(4)
1749     C endif
1750    
1751     c if (dPhy(np2).gt.dPhy(np4).and.dPhy(np4).gt.0. _d 0) then
1752     c print*,'QQQ dphy',np2,' > dphy',np4,i,j,k,Phy2(i,j,k),
1753     c & Ptr(i,j,k,bi,bj,iPhy+np4-1), dPhy(2), dPhy(4)
1754     endif
1755     c if (Ptr(i,j,k,bi,bj,iphy+np2-1).gt.Ptr(i,j,k,bi,bj,iPhy+np4-1)
1756     c & .and. Ptr(i,j,k,bi,bj,iPhy+np4-1).gt.0. _d 0) then
1757     c print*,'QQ phy',np2,' > ',np4,i,j,k,
1758     c & Ptr(i,j,k,bi,bj,iPhy+np2-1),
1759     c & Ptr(i,j,k,bi,bj,iPhy+np4-1), dPhy(2), dPhy(4)
1760     endif
1761    
1762     endif
1763     enddo ! np
1764     endif ! k
1765     endif
1766     #endif
1767    
1768     #ifdef ALLOW_DIAGNOSTICS
1769     COJ for diagnostics
1770     PParr(i,j,k) = PP
1771     Nfixarr(i,j,k) = Nfix
1772     c ANNA_TAVE
1773     #ifdef WAVES_DIAG_PCHL
1774     DO np = 1,npmax
1775     Pchlarr(i,j,k,np) = phychl(np)
1776     ENDDO
1777     #endif
1778     c ANNA end TAVE
1779     #ifdef DAR_DIAG_RSTAR
1780     DO np = 1,npmax
1781     Rstararr(i,j,k,np) = Rstarl(np)
1782     ENDDO
1783     #endif
1784     #ifdef ALLOW_DIAZ
1785     #ifdef DAR_DIAG_NFIXP
1786     DO np = 1,npmax
1787     NfixParr(i,j,k,np) = NfixPl(np)
1788     ENDDO
1789     #endif
1790     #endif
1791     #ifdef DAR_DIAG_CHL
1792     GeiderChlarr(i,j,k) = ChlGeiderlocal
1793     DoneyChlarr(i,j,k) = ChlDoneylocal
1794     CloernChlarr(i,j,k) = ChlCloernlocal
1795     IF (totphyC .NE. 0. _d 0) THEN
1796     GeiderChl2Carr(i,j,k) = ChlGeiderlocal/totphyC
1797     DoneyChl2Carr(i,j,k) = ChlDoneylocal/totphyC
1798     CloernChl2Carr(i,j,k) = ChlCloernlocal/totphyC
1799     ELSE
1800     GeiderChl2Carr(i,j,k) = 0. _d 0
1801     DoneyChl2Carr(i,j,k) = 0. _d 0
1802     CloernChl2Carr(i,j,k) = 0. _d 0
1803     ENDIF
1804     #endif
1805     COJ
1806     #endif /* ALLOW_DIAGNOSTICS */
1807    
1808     c total fixation (NOTE - STILL NEEDS GLOB SUM)
1809     tot_Nfix=tot_Nfix+
1810     & Nfix*rA(i,j,bi,bj)*rF(k)*hFacC(i,j,k,bi,bj)
1811    
1812     #ifdef ALLOW_TIMEAVE
1813     c save averages
1814     c Phygrow1ave(i,j,k,bi,bj)=Phygrow1ave(i,j,k,bi,bj)+
1815     c & mu1*py1*deltaTclock
1816     c & /float(nsubtime)
1817     c Phygrow2ave(i,j,k,bi,bj)=Phygrow2ave(i,j,k,bi,bj)+
1818     c & mu2*py2*deltaTclock
1819     c & /float(nsubtime)
1820     c Zoograzave(i,j,k,bi,bj)=Zoograzave(i,j,k,bi,bj)+
1821     c & (gampn1*graz1*zo +gampn2*graz2*zo)*
1822     c & deltaTclock/float(nsubtime)
1823     #ifdef GEIDER
1824     Chlave(i,j,k,bi,bj)=Chlave(i,j,k,bi,bj)+
1825     & Chl*dtplankton
1826     #endif
1827     PARave(i,j,k,bi,bj)=PARave(i,j,k,bi,bj)+
1828     & PARl*dtplankton
1829     PPave(i,j,k,bi,bj)=PPave(i,j,k,bi,bj)+
1830     & PP*dtplankton
1831     Nfixave(i,j,k,bi,bj)=Nfixave(i,j,k,bi,bj)+
1832     & Nfix*dtplankton
1833     Denitave(i,j,k,bi,bj)=Denitave(i,j,k,bi,bj)+
1834     & denit*dtplankton
1835     #ifdef WAVES_DIAG_PCHL
1836     do np=1,npmax
1837     Pchlave(i,j,k,bi,bj,np)=Pchlave(i,j,k,bi,bj,np)+
1838     & phychl(np)*dtplankton
1839     enddo
1840     #endif
1841     #ifdef DAR_DIAG_ACDOM
1842     c print*,'acdom',k,acdom_k(k,darwin_diag_acdom_ilam)
1843     aCDOMave(i,j,k,bi,bj)=aCDOMave(i,j,k,bi,bj)+
1844     & acdom_k(k,darwin_diag_acdom_ilam)*dtplankton
1845     #endif
1846     #ifdef DAR_DIAG_IRR
1847     do ilam = 1,tlam
1848     if (k.EQ.1) then
1849     Edave(i,j,k,bi,bj,ilam)=Edave(i,j,k,bi,bj,ilam)+
1850     & Edwsf(ilam)*dtplankton
1851     Esave(i,j,k,bi,bj,ilam)=Esave(i,j,k,bi,bj,ilam)+
1852     & Eswsf(ilam)*dtplankton
1853     Coj no Eu at surface (yet)
1854     else
1855     Edave(i,j,k,bi,bj,ilam)=Edave(i,j,k,bi,bj,ilam)+
1856     & Edz(ilam,k-1)*dtplankton
1857     Esave(i,j,k,bi,bj,ilam)=Esave(i,j,k,bi,bj,ilam)+
1858     & Esz(ilam,k-1)*dtplankton
1859     Euave(i,j,k,bi,bj,ilam)=Euave(i,j,k,bi,bj,ilam)+
1860     & Euz(ilam,k-1)*dtplankton
1861     endif
1862     Eutave(i,j,k,bi,bj,ilam)=Eutave(i,j,k,bi,bj,ilam)+
1863     & Eutop(ilam,k)*dtplankton
1864     enddo
1865     #endif
1866     #ifdef DAR_DIAG_ABSORP
1867     do ilam = 1,tlam
1868     aave(i,j,k,bi,bj,ilam)=aave(i,j,k,bi,bj,ilam)+
1869     & a_k(k,ilam)*dtplankton
1870     enddo
1871     #endif
1872     #ifdef DAR_DIAG_SCATTER
1873     do ilam = 1,tlam
1874     btave(i,j,k,bi,bj,ilam)=btave(i,j,k,bi,bj,ilam)+
1875     & bt_k(k,ilam)*dtplankton
1876     bbave(i,j,k,bi,bj,ilam)=bbave(i,j,k,bi,bj,ilam)+
1877     & bb_k(k,ilam)*dtplankton
1878     enddo
1879     #endif
1880     #ifdef DAR_DIAG_PART_SCATTER
1881     do ilam = 1,tlam
1882     apartave(i,j,k,bi,bj,ilam)=apartave(i,j,k,bi,bj,ilam)+
1883     & apart_k(k,ilam)*dtplankton
1884     btpartave(i,j,k,bi,bj,ilam)=btpartave(i,j,k,bi,bj,ilam)+
1885     & bpart_k(k,ilam)*dtplankton
1886     bbpartave(i,j,k,bi,bj,ilam)=bbpartave(i,j,k,bi,bj,ilam)+
1887     & bbpart_k(k,ilam)*dtplankton
1888     enddo
1889     #endif
1890     #ifdef DAR_DIAG_RSTAR
1891     do np=1,npmax
1892     Rstarave(i,j,k,bi,bj,np)=Rstarave(i,j,k,bi,bj,np)+
1893     & Rstarl(np)*dtplankton
1894     RNstarave(i,j,k,bi,bj,np)=RNstarave(i,j,k,bi,bj,np)+
1895     & RNstarl(np)*dtplankton
1896     enddo
1897     #endif
1898     #ifdef DAR_DIAG_DIVER
1899     Diver1ave(i,j,k,bi,bj)=Diver1ave(i,j,k,bi,bj)+
1900     & Diver1(i,j,k)*dtplankton
1901     Diver2ave(i,j,k,bi,bj)=Diver2ave(i,j,k,bi,bj)+
1902     & Diver2(i,j,k)*dtplankton
1903     Diver3ave(i,j,k,bi,bj)=Diver3ave(i,j,k,bi,bj)+
1904     & Diver3(i,j,k)*dtplankton
1905     Diver4ave(i,j,k,bi,bj)=Diver4ave(i,j,k,bi,bj)+
1906     & Diver4(i,j,k)*dtplankton
1907     #endif
1908     #ifdef DAR_DIAG_GROW
1909     do np=1,npmax
1910     Growave(i,j,k,bi,bj,np)=Growave(i,j,k,bi,bj,np)+
1911     & Growl(np)*dtplankton
1912     Growsqave(i,j,k,bi,bj,np)=Growsqave(i,j,k,bi,bj,np)+
1913     & Growsql(np)*dtplankton
1914     enddo
1915     #endif
1916    
1917     #ifdef ALLOW_DIAZ
1918     #ifdef DAR_DIAG_NFIXP
1919     do np=1,npmax
1920     NfixPave(i,j,k,bi,bj,np)=NfixPave(i,j,k,bi,bj,np)+
1921     & NfixPl(np)*dtplankton
1922     enddo
1923     #endif
1924     #endif
1925     #endif
1926    
1927     #ifdef ALLOW_CARBON
1928     if (k.eq.1) then
1929     SURave(i,j,bi,bj) =SURave(i,j,bi,bj)+
1930     & flxCO2(i,j)*dtplankton
1931     SURCave(i,j,bi,bj) =SURCave(i,j,bi,bj)+
1932     & FluxCO2(i,j,bi,bj)*dtplankton
1933     SUROave(i,j,bi,bj) =SUROave(i,j,bi,bj)+
1934     & flxO2(i,j)*dtplankton
1935     pCO2ave(i,j,bi,bj) =pCO2ave(i,j,bi,bj)+
1936     & pCO2(i,j,bi,bj)*dtplankton
1937     pHave(i,j,bi,bj) =pHave(i,j,bi,bj)+
1938     & pH(i,j,bi,bj)*dtplankton
1939     endif
1940     #endif
1941     endif
1942     c end if hFac>0
1943    
1944     enddo ! k
1945     c end layer loop
1946     c
1947     ENDDO ! i
1948     ENDDO ! j
1949    
1950     #ifdef CO2_FLUX_BUDGET
1951     C find k index from maximum GGL90 mixing length,
1952     C use this as "mixing layer depth"
1953     DO k=1,Nr
1954     DO j=jmin,jmax
1955     DO i=imin,imax
1956     if(mixingLength(i,j,k,bi,bj) .EQ.
1957     & MAXVAL(mixingLength(i,j,1:Nr,bi,bj))) then
1958     mixingDepthKLev(i,j) = MIN(k,Nr)
1959     mixingDepth(i,j) = ABS(rF(MIN(k,Nr)))
1960     endif
1961     ENDDO
1962     ENDDO
1963     ENDDO
1964    
1965     C vertically-integrate relevant biological DIC tendency terms
1966     C mmol C m^-3 s^-1
1967     DO k=1,Nr
1968     DO j=jmin,jmax
1969     DO i=imin,imax
1970    
1971     if(k .LE. mixingDepthKLev(i,j)) then
1972     deltaDic_bio(i,j) = deltaDic_bio(i,j) +
1973     & (-(budgetConsumpDIC(i,j,k,bi,bj)) +
1974     & -(budgetConsumpDIC_PIC(i,j,k,bi,bj)) +
1975     & budgetDOCRemin(i,j,k,bi,bj) +
1976     & budgetPReminC(i,j,k,bi,bj) +
1977     & disscPIC(i,j,k,bi,bj))
1978     endif
1979    
1980     ENDDO
1981     ENDDO
1982     ENDDO
1983    
1984     C compute CO2 flux budget terms (mmol C m^-2 s^-1)
1985     DO j=jmin,jmax
1986     DO i=imin,imax
1987    
1988     dCO2Flux_theta(i,j) = deltaDic_theta(i,j) *
1989     & mixingDepth(i,j) / deltaT *
1990     & (1. _d 0 - FIce(i,j,bi,bj))
1991    
1992     dCO2Flux_salt(i,j) = deltaDic_salt(i,j) *
1993     & mixingDepth(i,j) / deltaT *
1994     & (1. _d 0 - FIce(i,j,bi,bj))
1995    
1996     dCO2Flux_alk(i,j) = deltaDic_alk(i,j) *
1997     & mixingDepth(i,j) / deltaT *
1998     & (1. _d 0 - FIce(i,j,bi,bj))
1999    
2000     dCO2Flux_apCO2(i,j) = deltaDic_apCO2(i,j) *
2001     & mixingDepth(i,j) / deltaT *
2002     & (1. _d 0 - FIce(i,j,bi,bj))
2003    
2004     dCO2Flux_residual(i,j) = deltaDic_residual(i,j) *
2005     & mixingDepth(i,j) / deltaT *
2006     & (1. _d 0 - FIce(i,j,bi,bj))
2007    
2008     dCO2Flux_bio(i,j) = deltaDic_bio(i,j) *
2009     & mixingDepth(i,j) *
2010     & (1. _d 0 - FIce(i,j,bi,bj))
2011    
2012     dCO2Flux_circ(i,j) =
2013     & dCO2Flux_residual(i,j) -
2014     & dCO2Flux_bio(i,j)
2015    
2016     ENDDO
2017     ENDDO
2018    
2019     if (budgetTStep1 .EQ. 0) then
2020     budgetTStep1 = 1;
2021     endif
2022    
2023     #endif /* CO2_FLUX_BUDGET */
2024    
2025     #ifdef ALLOW_PAR_DAY
2026     C 1 <-> 2
2027     PARiaccum = 3 - PARiprev
2028    
2029     DO k=1,nR
2030     DO j=1,sNy
2031     DO i=1,sNx
2032     PARday(i,j,k,bi,bj,PARiaccum) =
2033     & PARday(i,j,k,bi,bj,PARiaccum) + PAR(i,j,k)
2034     ENDDO
2035     ENDDO
2036     ENDDO
2037    
2038     phase = 0. _d 0
2039     itistime = DIFF_PHASE_MULTIPLE( phase, darwin_PARavPeriod,
2040     & newtime, dtsubtime)
2041    
2042     IF ( itistime ) THEN
2043     C compute average
2044     nav = darwin_PARnav
2045     IF (newtime - baseTime .LT. darwin_PARavPeriod) THEN
2046     C incomplete period at beginning of run
2047     nav = NINT((newtime-baseTime)/dtsubtime)
2048     ENDIF
2049     DO k=1,nR
2050     DO j=1,sNy
2051     DO i=1,sNx
2052     PARday(i,j,k,bi,bj,PARiaccum) =
2053     & PARday(i,j,k,bi,bj,PARiaccum) / nav
2054     ENDDO
2055     ENDDO
2056     ENDDO
2057     C reset the other slot for averaging
2058     DO k=1,nR
2059     DO j=1,sNy
2060     DO i=1,sNx
2061     PARday(i,j,k,bi,bj,PARiprev) = 0. _d 0
2062     ENDDO
2063     ENDDO
2064     ENDDO
2065     ENDIF
2066     C itistime
2067     #endif
2068    
2069     COJ fill diagnostics
2070     #ifdef ALLOW_DIAGNOSTICS
2071     IF ( useDiagnostics ) THEN
2072     diagname = ' '
2073     WRITE(diagname,'(A8)') 'PAR '
2074     CALL DIAGNOSTICS_FILL( PAR(1-Olx,1-Oly,1), diagname,
2075     & 0,Nr,2,bi,bj,myThid )
2076     WRITE(diagname,'(A8)') 'PP '
2077     CALL DIAGNOSTICS_FILL( PParr(1-Olx,1-Oly,1), diagname,
2078     & 0,Nr,2,bi,bj,myThid )
2079     WRITE(diagname,'(A8)') 'Nfix '
2080     CALL DIAGNOSTICS_FILL( Nfixarr(1-Olx,1-Oly,1), diagname,
2081     & 0,Nr,2,bi,bj,myThid )
2082     c ANNA_TAVE
2083     #ifdef WAVES_DIAG_PCHL
2084     DO np=1,MIN(99,npmax)
2085     WRITE(diagname,'(A5,I2.2,A1)') 'Pchl',np,' '
2086     CALL DIAGNOSTICS_FILL( Pchlarr(1-Olx,1-Oly,1,np), diagname,
2087     & 0,Nr,2,bi,bj,myThid )
2088     ENDDO
2089     #endif
2090     c ANNA end TAVE
2091     #ifdef DAR_DIAG_RSTAR
2092     DO np=1,MIN(99,npmax)
2093     WRITE(diagname,'(A5,I2.2,A1)') 'Rstar',np,' '
2094     CALL DIAGNOSTICS_FILL( Rstararr(1-Olx,1-Oly,1,np), diagname,
2095     & 0,Nr,2,bi,bj,myThid )
2096     ENDDO
2097     #endif
2098     #ifdef DAR_DIAG_DIVER
2099     WRITE(diagname,'(A8)') 'Diver1 '
2100     CALL DIAGNOSTICS_FILL( Diver1(1-Olx,1-Oly,1), diagname,
2101     & 0,Nr,2,bi,bj,myThid )
2102     WRITE(diagname,'(A8)') 'Diver2 '
2103     CALL DIAGNOSTICS_FILL( Diver2(1-Olx,1-Oly,1), diagname,
2104     & 0,Nr,2,bi,bj,myThid )
2105     WRITE(diagname,'(A8)') 'Diver3 '
2106     CALL DIAGNOSTICS_FILL( Diver3(1-Olx,1-Oly,1), diagname,
2107     & 0,Nr,2,bi,bj,myThid )
2108     WRITE(diagname,'(A8)') 'Diver4 '
2109     CALL DIAGNOSTICS_FILL( Diver4(1-Olx,1-Oly,1), diagname,
2110     & 0,Nr,2,bi,bj,myThid )
2111     #endif
2112     #ifdef ALLOW_DIAZ
2113     #ifdef DAR_DIAG_NFIXP
2114     DO np=1,MIN(99,npmax)
2115     WRITE(diagname,'(A5,I2.2,A1)') 'NfixP',np,' '
2116     CALL DIAGNOSTICS_FILL( NfixParr(1-Olx,1-Oly,1,np), diagname,
2117     & 0,Nr,2,bi,bj,myThid )
2118     ENDDO
2119     #endif
2120     #endif
2121     #ifdef DAR_DIAG_CHL
2122     CALL DIAGNOSTICS_FILL( GeiderChlarr(1-Olx,1-Oly,1), 'ChlGeide',
2123     & 0,Nr,2,bi,bj,myThid )
2124     CALL DIAGNOSTICS_FILL( GeiderChl2Carr(1-Olx,1-Oly,1),'Chl2CGei',
2125     & 0,Nr,2,bi,bj,myThid )
2126     CALL DIAGNOSTICS_FILL( DoneyChlarr(1-Olx,1-Oly,1), 'ChlDoney',
2127     & 0,Nr,2,bi,bj,myThid )
2128     CALL DIAGNOSTICS_FILL( DoneyChl2Carr(1-Olx,1-Oly,1), 'Chl2CDon',
2129     & 0,Nr,2,bi,bj,myThid )
2130     CALL DIAGNOSTICS_FILL( CloernChlarr(1-Olx,1-Oly,1), 'ChlCloer',
2131     & 0,Nr,2,bi,bj,myThid )
2132     CALL DIAGNOSTICS_FILL( CloernChl2Carr(1-Olx,1-Oly,1),'Chl2CClo',
2133     & 0,Nr,2,bi,bj,myThid )
2134     #endif
2135     #ifdef ALLOW_CARBON
2136     CALL DIAGNOSTICS_FILL( flxCO2(1-Olx,1-Oly), 'DICTFLX ',
2137     & 0,1,2,bi,bj,myThid )
2138     CALL DIAGNOSTICS_FILL( FluxCO2(1-Olx,1-Oly,bi,bj), 'DICCFLX ',
2139     & 0,1,2,bi,bj,myThid )
2140     CALL DIAGNOSTICS_FILL( flxO2(1-Olx,1-Oly), 'DICOFLX ',
2141     & 0,1,2,bi,bj,myThid )
2142     CALL DIAGNOSTICS_FILL( fugf(1-Olx,1-Oly,bi,bj), 'DICFGCO2',
2143     & 0,1,2,bi,bj,myThid )
2144     CALL DIAGNOSTICS_FILL( pCO2(1-Olx,1-Oly,bi,bj), 'DICPCO2 ',
2145     & 0,1,2,bi,bj,myThid )
2146     CALL DIAGNOSTICS_FILL( pH(1-Olx,1-Oly,bi,bj), 'DICPHAV ',
2147     & 0,1,2,bi,bj,myThid )
2148     CALL DIAGNOSTICS_FILL(KspTP(1-Olx,1-Oly,bi,bj), 'KSPTP ',
2149     & 0,1,2,bi,bj,myThid)
2150     CALL DIAGNOSTICS_FILL(calcium(1-Olx,1-Oly,bi,bj), 'CALCIUM ',
2151     & 0,1,2,bi,bj,myThid)
2152     CALL DIAGNOSTICS_FILL(omegaC(1-Olx,1-Oly,1,bi,bj), 'OMEGAC ',
2153     & 0,Nr,2,bi,bj,myThid)
2154     CALL DIAGNOSTICS_FILL(disscPIC(1-Olx,1-Oly,1,bi,bj), 'DISSC ',
2155     & 0,Nr,2,bi,bj,myThid)
2156     #ifdef ALLOW_SED_DISS_FLUX
2157     CALL DIAGNOSTICS_FILL(DICSedFlux(1-Olx,1-Oly,bi,bj), 'DICSFLX ',
2158     & 0,1,2,bi,bj,myThid)
2159     CALL DIAGNOSTICS_FILL(ALKSedFlux(1-Olx,1-Oly,bi,bj), 'ALKSFLX ',
2160     & 0,1,2,bi,bj,myThid)
2161     CALL DIAGNOSTICS_FILL(RFlux(1-Olx,1-Oly,bi,bj), 'RFLUX ',
2162     & 0,1,2,bi,bj,myThid)
2163     CALL DIAGNOSTICS_FILL(CO3Sw(1-Olx,1-Oly,bi,bj), 'CO3SW ',
2164     & 0,1,2,bi,bj,myThid)
2165     CALL DIAGNOSTICS_FILL(CO3Sed(1-Olx,1-Oly,bi,bj), 'CO3SED ',
2166     & 0,1,2,bi,bj,myThid)
2167     #endif /* ALLOW_SED_DISS_FLUX */
2168     #ifdef CO2_FLUX_BUDGET
2169     CALL DIAGNOSTICS_FILL(dCO2Flux(1-Olx,1-Oly,bi,bj),
2170     & 'DCO2FLX ',0,1,2,bi,bj,myThid )
2171     CALL DIAGNOSTICS_FILL(dCO2Flux_theta(1-Olx,1-Oly),
2172     & 'DCO2FLXT',0,1,2,bi,bj,myThid)
2173     CALL DIAGNOSTICS_FILL(dCO2Flux_salt(1-Olx,1-Oly),
2174     & 'DCO2FLXS',0,1,2,bi,bj,myThid)
2175     CALL DIAGNOSTICS_FILL(dCO2Flux_alk(1-Olx,1-Oly),
2176     & 'DCO2FLXA',0,1,2,bi,bj,myThid)
2177     CALL DIAGNOSTICS_FILL(dCO2Flux_apCO2(1-Olx,1-Oly),
2178     & 'DCO2FLXC',0,1,2,bi,bj,myThid)
2179     CALL DIAGNOSTICS_FILL(dCO2Flux_residual(1-Olx,1-Oly),
2180     & 'DCO2FLXR',0,1,2,bi,bj,myThid)
2181     CALL DIAGNOSTICS_FILL(dCO2Flux_bio(1-Olx,1-Oly),
2182     & 'DCO2FLXB',0,1,2,bi,bj,myThid)
2183     CALL DIAGNOSTICS_FILL(dCO2Flux_circ(1-Olx,1-Oly),
2184     & 'DCO2FLXP',0,1,2,bi,bj,myThid)
2185     CALL DIAGNOSTICS_FILL(mixingDepth(1-Olx,1-Oly),
2186     & 'BCO2MIXD',0,1,2,bi,bj,myThid)
2187     #endif /* CO2_FLUX_BUDGET */
2188     #endif /* ALLOW_CARBON */
2189     ENDIF
2190     #endif /* ALLOW_DIAGNOSTICS */
2191     COJ
2192    
2193     c determine iron partitioning - solve for free iron
2194     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
2195     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
2196     & myIter, mythid)
2197     c
2198     #ifdef ALLOW_TIMEAVE
2199     c save averages
2200     dar_timeave(bi,bj) = dar_timeave(bi,bj) + dtplankton
2201     #ifdef ALLOW_CARBON
2202     dic_timeave(bi,bj) = dic_timeave(bi,bj) + dtplankton
2203     #endif
2204     #endif
2205     c
2206     c -----------------------------------------------------
2207     ENDDO ! it
2208     c -----------------------------------------------------
2209     c end of bio-chemical time loop
2210     c
2211     RETURN
2212     END
2213     #endif /*DARWIN*/
2214     #endif /*ALLOW_PTRACERS*/
2215    
2216     C============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22