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