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