/[MITgcm]/MITgcm/pkg/dic/dic_biotic_forcing.F
ViewVC logotype

Contents of /MITgcm/pkg/dic/dic_biotic_forcing.F

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


Revision 1.32 - (show annotations) (download)
Sat Jan 16 21:57:12 2016 UTC (8 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, HEAD
Changes since 1.31: +9 -1 lines
- use the updated seaice fraction (from seaice pkgs) for DIC forcing
  instead of the previous time-step value ; done in new S/R DIC_FIELDS_UPDATE
  called from the top of dic_biotic_forcing.F (instead of in DIC_FIELDS_LOAD);
- move also update from ATM-OCN coupler in the same new S/R (previously done
  in ocn_apply_import.F).

1 C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_biotic_forcing.F,v 1.31 2014/12/08 22:59:04 jmc Exp $
2 C $Name: $
3
4 #include "DIC_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: DIC_BIOTIC_FORCING
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE DIC_BIOTIC_FORCING(
11 U PTR_DIC, PTR_ALK, PTR_PO4, PTR_DOP,
12 #ifdef ALLOW_O2
13 U PTR_O2,
14 #endif
15 #ifdef ALLOW_FE
16 U PTR_FE,
17 #endif
18 I bi, bj, iMin, iMax, jMin, jMax,
19 I myIter, myTime, myThid )
20
21 C !DESCRIPTION:
22 C updates all the tracers for the effects of air-sea exchange, biological
23 c activity and remineralization
24
25 C !USES: ===============================================================
26 IMPLICIT NONE
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31 #include "DYNVARS.h"
32 #include "DIC_VARS.h"
33 #include "PTRACERS_SIZE.h"
34 #include "PTRACERS_PARAMS.h"
35
36 C !INPUT/OUTPUT PARAMETERS: ===================================================
37 C PTR_DIC :: dissolced inorganic carbon
38 C PTR_ALK :: alkalinity
39 C PTR_PO4 :: phosphate
40 c PTR_DOP :: dissolve organic phosphurous
41 c PTR_O2 :: oxygen
42 C PTR_FE :: iron
43 c bi, bj :: current tile indices
44 C myIter :: current timestep
45 C myTime :: current time
46 C myThid :: thread number
47 _RL PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48 _RL PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49 _RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
50 _RL PTR_DOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51 #ifdef ALLOW_O2
52 _RL PTR_O2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
53 #endif
54 #ifdef ALLOW_FE
55 _RL PTR_FE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
56 #endif
57 INTEGER bi, bj, iMin, iMax, jMin, jMax
58 INTEGER myIter
59 _RL myTime
60 INTEGER myThid
61
62 #ifdef ALLOW_PTRACERS
63 #ifdef DIC_BIOTIC
64
65 C !LOCAL VARIABLES: ====================================================
66 C i,j,k :: loop indices
67 C G* :: tendency term for the tracers
68 C SURA :: tendency of alkalinity due to freshwater
69 C SURC :: tendency of DIC due to air-sea exchange
70 C and virtual flux
71 C SURO :: tendency of O2 due to air-sea exchange
72 C GPO4 :: tendency of PO4 due to biological productivity,
73 C exchange with DOP pool and reminerization
74 C CAR :: carbonate changes due to biological
75 C productivity and remineralization
76 C BIOac :: biological productivity
77 C RDOP :: DOP sink due to remineralization
78 C pflux :: changes to PO4 due to flux and remineralization
79 C CAR_S :: carbonate sink
80 C cflux :: carbonate changes due to flux and remineralization
81 C freefe :: iron not bound to ligand
82 _RL GDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83 _RL GALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84 _RL GPO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85 _RL GDOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86 _RL SURA(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL SURC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL SURO(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL CAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90 _RL BIOac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91 _RL RDOP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92 _RL pflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93 _RL exportflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94 _RL CAR_S(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95 _RL cflux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96 #ifdef ALLOW_O2
97 _RL GO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98 #endif
99 #ifdef ALLOW_FE
100 _RL GFE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101 _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102 #endif
103 INTEGER i,j,k
104 #ifdef CAR_DISS
105 INTEGER nCALCITEstep
106 #endif
107 #ifdef ALLOW_FE
108 # ifdef SEDFE
109 INTEGER kBottom
110 # endif
111 #endif
112 CEOP
113
114 #ifdef ALLOW_DEBUG
115 IF (debugMode) CALL DEBUG_ENTER('DIC_BIOTIC_FORCING',myThid)
116 #endif
117
118 IF ( useThSIce .OR. useSEAICE .OR. useCoupler ) THEN
119 #ifdef ALLOW_DEBUG
120 IF (debugMode) CALL DEBUG_CALL('DIC_FIELDS_UPDATE',myThid)
121 #endif
122 CALL DIC_FIELDS_UPDATE(
123 I bi, bj, myTime, myIter, myThid )
124 ENDIF
125
126 DO k=1,Nr
127 DO j=1-OLy,sNy+OLy
128 DO i=1-OLx,sNx+OLx
129 RDOP(i,j,k) =0. _d 0
130 GDIC(i,j,k) =0. _d 0
131 GALK(i,j,k) =0. _d 0
132 GPO4(i,j,k) =0. _d 0
133 GDOP(i,j,k) =0. _d 0
134 CAR(i,j,k) =0. _d 0
135 BIOac(i,j,k) =0. _d 0
136 pflux(i,j,k) =0. _d 0
137 exportflux(i,j,k)=0. _d 0
138 cflux(i,j,k) =0. _d 0
139 CAR_S(i,j,k) =0. _d 0
140 #ifdef ALLOW_O2
141 GO2(i,j,k) =0. _d 0
142 #endif
143 #ifdef ALLOW_FE
144 GFE(i,j,k) =0. _d 0
145 C no longer needed after adding full initialisation of freefe in S/R FE_CHEM
146 c freefe(i,j,k) =0. _d 0
147 #endif
148 ENDDO
149 ENDDO
150 ENDDO
151 DO j=1-OLy,sNy+OLy
152 DO i=1-OLx,sNx+OLx
153 SURA(i,j) =0. _d 0
154 SURC(i,j) =0. _d 0
155 SURO(i,j) =0. _d 0
156 ENDDO
157 ENDDO
158
159 C carbon air-sea interaction
160 #ifdef ALLOW_DEBUG
161 IF (debugMode) CALL DEBUG_CALL('DIC_SURFFORCING',myThid)
162 #endif
163 CALL DIC_SURFFORCING(
164 I PTR_DIC, PTR_ALK, PTR_PO4,
165 O SURC,
166 I bi, bj, iMin, iMax, jMin, jMax,
167 I myIter, myTime, myThid )
168
169 C alkalinity air-sea interaction
170 #ifdef ALLOW_DEBUG
171 IF (debugMode) CALL DEBUG_CALL('ALK_SURFFORCING',myThid)
172 #endif
173 CALL ALK_SURFFORCING(
174 I PTR_ALK,
175 O SURA,
176 I bi, bj, iMin, iMax, jMin, jMax,
177 I myIter, myTime, myThid )
178
179 #ifdef ALLOW_O2
180 C oxygen air-sea interaction
181 #ifdef ALLOW_DEBUG
182 IF (debugMode) CALL DEBUG_CALL('O2_SURFFORCING',myThid)
183 #endif
184 CALL O2_SURFFORCING(
185 I PTR_O2,
186 O SURO,
187 I bi, bj, iMin, iMax, jMin, jMax,
188 I myIter, myTime, myThid )
189 #endif
190
191 #ifdef ALLOW_FE
192 C find free iron
193 #ifdef ALLOW_DEBUG
194 IF (debugMode) CALL DEBUG_CALL('FE_CHEM',myThid)
195 #endif
196 CALL FE_CHEM( bi, bj, iMin, iMax, jMin, jMax,
197 U PTR_FE,
198 O freefe,
199 I myIter, myThid )
200 #endif
201
202 C biological activity
203 #ifdef ALLOW_DEBUG
204 IF (debugMode) CALL DEBUG_CALL('BIO_EXPORT',myThid)
205 #endif
206 CALL BIO_EXPORT(
207 I PTR_PO4,
208 #ifdef ALLOW_FE
209 I PTR_FE,
210 #endif
211 O BIOac,
212 I bi, bj, iMin, iMax, jMin, jMax,
213 I myIter, myTime, myThid )
214
215 C flux of po4 from layers with biological activity
216 #ifdef ALLOW_DEBUG
217 IF (debugMode) CALL DEBUG_CALL('PHOS_FLUX',myThid)
218 #endif
219 CALL PHOS_FLUX(
220 I BIOac,
221 U pflux, exportflux,
222 I bi, bj, iMin, iMax, jMin, jMax,
223 I myIter, myTime, myThid )
224
225 C- Carbonate sink
226 DO k=1,Nr
227 DO j=jMin,jMax
228 DO i=iMin,iMax
229 CAR_S(i,j,k)=BIOac(i,j,k)*R_CP*rain_ratio(i,j,bi,bj)*
230 & (1. _d 0-DOPfraction)
231 ENDDO
232 ENDDO
233 ENDDO
234
235 C carbonate
236 #ifdef CAR_DISS
237 C dissolution only below saturation horizon
238 C code following method by Karsten Friis
239 nCALCITEstep = 3600
240 IF(myIter .lt. (nIter0+5) .or.
241 & mod(myIter,nCALCITEstep) .eq. 0)THEN
242 #ifdef ALLOW_DEBUG
243 IF (debugMode) CALL DEBUG_CALL('CALCITE_SATURATION',myThid)
244 #endif
245 CALL CALCITE_SATURATION(
246 I PTR_DIC, PTR_ALK, PTR_PO4,
247 I bi, bj, iMin, iMax, jMin, jMax,
248 I myIter, myTime, myThid )
249 ENDIF
250
251 #ifdef ALLOW_DEBUG
252 IF (debugMode) CALL DEBUG_CALL('CAR_FLUX_OMEGA_TOP',myThid)
253 #endif
254 CALL CAR_FLUX_OMEGA_TOP(
255 I BIOac,
256 O cflux,
257 I bi, bj, iMin, iMax, jMin, jMax,
258 I myIter, myTime, myThid )
259 #else
260 C old OCMIP way
261 #ifdef ALLOW_DEBUG
262 IF (debugMode) CALL DEBUG_CALL('CAR_FLUX',myThid)
263 #endif
264 CALL CAR_FLUX(
265 I CAR_S,
266 U cflux,
267 I bi, bj, iMin, iMax, jMin, jMax,
268 I myIter, myTime, myThid )
269 #endif
270
271 C add all tendencies for PO4, DOP, ALK, DIC
272 DO k=1,Nr
273 DO j=jMin,jMax
274 DO i=iMin,iMax
275 #ifdef DIC_NO_NEG
276 RDOP(i,j,k)= MAX(maskC(i,j,k,bi,bj)*KDOPRemin*PTR_DOP(i,j,k)
277 & ,0. _d 0)
278 #else
279 RDOP(i,j,k)= maskC(i,j,k,bi,bj)*KDOPRemin*PTR_DOP(i,j,k)
280 #endif
281 GPO4(i,j,k)=-BIOac(i,j,k)+pflux(i,j,k) + RDOP(i,j,k)
282
283 car(i,j,k) = cflux(i,j,k) - CAR_S(i,j,k)
284
285 GDOP(i,j,k)=+BIOac(i,j,k)*DOPfraction - RDOP(i,j,k)
286
287 GALK(i,j,k)=+2. _d 0 *car(i,j,k)-R_NP*GPO4(i,j,k)
288
289 GDIC(i,j,k)=car(i,j,k)+R_CP*GPO4(i,j,k)
290
291 #ifdef ALLOW_O2
292 if (PTR_O2(i,j,k).GT.O2crit) then
293 GO2(i,j,k)= R_OP*GPO4(i,j,k)
294 else
295 GO2(i,j,k)= 0. _d 0
296 endif
297 #endif
298 #ifdef ALLOW_FE
299 GFE(i,j,k) = R_FeP*GPO4(i,j,k)
300 & -Kscav*freefe(i,j,k)
301 #endif
302 ENDDO
303 ENDDO
304 ENDDO
305
306 DO j=jMin,jMax
307 DO i=iMin,iMax
308 GALK(i,j,1)=GALK(i,j,1)+SURA(i,j)
309 GDIC(i,j,1)=GDIC(i,j,1)+SURC(i,j)
310 #ifdef ALLOW_O2
311 GO2(i,j,1) =GO2(i,j,1)+SURO(i,j)
312 #endif
313 #ifdef ALLOW_FE
314 GFE(i,j,1)=GFE(i,j,1)+alpfe*
315 & InputFe(i,j,bi,bj)*recip_drF(1)
316 & *recip_hFacC(i,j,1,bi,bj)
317 # ifdef SEDFE
318 C include iron sediment source using the flux of po4 into bottom layer
319 kBottom = MAX(kLowC(i,j,bi,bj),1)
320 GFE(i,j,kBottom)=GFE(i,j,kBottom)
321 & +( fesedflux_pcm*pflux(i,j,kBottom) + FeIntSec )
322 & *recip_drF(kBottom)*recip_hFacC(i,j,kBottom,bi,bj)
323 # endif
324 #endif
325 ENDDO
326 ENDDO
327
328 IF ( useOBCS ) THEN
329 DO k=1,Nr
330 DO j=jMin,jMax
331 DO i=iMin,iMax
332 GDIC(i,j,k) = GDIC(i,j,k)*maskInC(i,j,bi,bj)
333 GALK(i,j,k) = GALK(i,j,k)*maskInC(i,j,bi,bj)
334 GPO4(i,j,k) = GPO4(i,j,k)*maskInC(i,j,bi,bj)
335 GDOP(i,j,k) = GDOP(i,j,k)*maskInC(i,j,bi,bj)
336 #ifdef ALLOW_O2
337 GO2(i,j,k) = GO2(i,j,k)*maskInC(i,j,bi,bj)
338 #endif
339 #ifdef ALLOW_FE
340 GFE(i,j,k) = GFE(i,j,k)*maskInC(i,j,bi,bj)
341 #endif
342 ENDDO
343 ENDDO
344 ENDDO
345 ENDIF
346
347 C update
348 DO k=1,Nr
349 DO j=jMin,jMax
350 DO i=iMin,iMax
351 PTR_DIC(i,j,k)=
352 & PTR_DIC(i,j,k)+GDIC(i,j,k)*PTRACERS_dTLev(k)
353 PTR_ALK(i,j,k)=
354 & PTR_ALK(i,j,k)+GALK(i,j,k)*PTRACERS_dTLev(k)
355 PTR_PO4(i,j,k)=
356 & PTR_PO4(i,j,k)+GPO4(i,j,k)*PTRACERS_dTLev(k)
357 PTR_DOP(i,j,k)=
358 & PTR_DOP(i,j,k)+GDOP(i,j,k)*PTRACERS_dTLev(k)
359 #ifdef ALLOW_O2
360 PTR_O2(i,j,k)=
361 & PTR_O2(i,j,k)+GO2(i,j,k)*PTRACERS_dTLev(k)
362 #endif
363 #ifdef ALLOW_FE
364 PTR_FE(i,j,k)=
365 & PTR_FE(i,j,k)+GFE(i,j,k)*PTRACERS_dTLev(k)
366 #endif
367 ENDDO
368 ENDDO
369 ENDDO
370
371 #ifdef ALLOW_FE
372 #ifdef MINFE
373 c find free iron and get rid of insoluble part
374 #ifdef ALLOW_DEBUG
375 IF (debugMode) CALL DEBUG_CALL('FE_CHEM',myThid)
376 #endif
377 CALL FE_CHEM( bi, bj, iMin, iMax, jMin, jMax,
378 U PTR_FE,
379 O freefe,
380 I myIter, myThid )
381 #endif
382 #endif
383
384 #ifdef ALLOW_TIMEAVE
385 C save averages
386 IF ( PTRACERS_taveFreq.GT.0. ) THEN
387 DO k=1,Nr
388 DO j=jMin,jMax
389 DO i=iMin,iMax
390 BIOave(i,j,k,bi,bj) =BIOave(i,j,k,bi,bj)+
391 & BIOac(i,j,k)*deltaTClock
392 CARave(i,j,k,bi,bj) =CARave(i,j,k,bi,bj)+
393 & CAR(i,j,k)*deltaTClock
394 OmegaCave(i,j,k,bi,bj)=OmegaCave(i,j,k,bi,bj)+
395 & OmegaC(i,j,k,bi,bj)*deltaTClock
396 pfluxave(i,j,k,bi,bj) =pfluxave(i,j,k,bi,bj) +
397 & pflux(i,j,k)*deltaTClock
398 epfluxave(i,j,k,bi,bj)=epfluxave(i,j,k,bi,bj) +
399 & exportflux(i,j,k)*deltaTClock
400 cfluxave(i,j,k,bi,bj) =cfluxave(i,j,k,bi,bj) +
401 & cflux(i,j,k)*deltaTClock
402 ENDDO
403 ENDDO
404 ENDDO
405 DO j=jMin,jMax
406 DO i=iMin,iMax
407 SURave(i,j,bi,bj) =SURave(i,j,bi,bj)+
408 & SURC(i,j)*deltaTClock
409 #ifdef ALLOW_O2
410 SUROave(i,j,bi,bj) =SUROave(i,j,bi,bj)+
411 & SURO(i,j)*deltaTClock
412 #endif
413 pCO2ave(i,j,bi,bj) =pCO2ave(i,j,bi,bj)+
414 & pCO2(i,j,bi,bj)*deltaTClock
415 pHave(i,j,bi,bj) =pHave(i,j,bi,bj)+
416 & pH(i,j,bi,bj)*deltaTClock
417 fluxCO2ave(i,j,bi,bj)=fluxCO2ave(i,j,bi,bj)+
418 & fluxCO2(i,j,bi,bj)*deltaTClock
419 ENDDO
420 ENDDO
421 DIC_timeAve(bi,bj) = DIC_timeAve(bi,bj)+deltaTClock
422 ENDIF
423 #endif /* ALLOW_TIMEAVE*/
424
425 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
426
427 #ifdef ALLOW_DIAGNOSTICS
428 IF ( useDiagnostics ) THEN
429 CALL DIAGNOSTICS_FILL(BIOac ,'DICBIOA ',0,Nr,2,bi,bj,myThid)
430 CALL DIAGNOSTICS_FILL(CAR ,'DICCARB ',0,Nr,2,bi,bj,myThid)
431 CALL DIAGNOSTICS_FILL(pCO2 ,'DICPCO2 ',0,1 ,1,bi,bj,myThid)
432 CALL DIAGNOSTICS_FILL(fluxCO2,'DICCFLX ',0,1 ,1,bi,bj,myThid)
433 CALL DIAGNOSTICS_FILL(pH ,'DICPHAV ',0,1 ,1,bi,bj,myThid)
434 CALL DIAGNOSTICS_FILL(SURC ,'DICTFLX ',0,1 ,2,bi,bj,myThid)
435 #ifdef ALLOW_O2
436 CALL DIAGNOSTICS_FILL(SURO ,'DICOFLX ',0,1 ,2,bi,bj,myThid)
437 #endif
438 ENDIF
439 #endif /* ALLOW_DIAGNOSTICS */
440
441 #ifdef ALLOW_DEBUG
442 IF (debugMode) CALL DEBUG_LEAVE('DIC_BIOTIC_FORCING',myThid)
443 #endif
444
445 #endif /* DIC_BIOTIC */
446 #endif /* ALLOW_PTRACERS */
447
448 RETURN
449 END

  ViewVC Help
Powered by ViewVC 1.1.22