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