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