/[MITgcm]/MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F
ViewVC logotype

Contents of /MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F

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


Revision 1.4 - (show annotations) (download)
Tue May 19 14:32:43 2015 UTC (10 years, 2 months ago) by benw
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt65w_20160512, ctrb_darwin2_ckpt66g_20170424, ctrb_darwin2_ckpt66k_20171025, ctrb_darwin2_ckpt66n_20180118, ctrb_darwin2_ckpt65v_20160409, ctrb_darwin2_ckpt65s_20160114, ctrb_darwin2_ckpt66d_20170214, ctrb_darwin2_ckpt65m_20150615, ctrb_darwin2_ckpt65q_20151118, ctrb_darwin2_ckpt65o_20150914, ctrb_darwin2_ckpt65p_20151023, ctrb_darwin2_ckpt65z_20160929, ctrb_darwin2_ckpt65n_20150729, ctrb_darwin2_ckpt66h_20170602, ctrb_darwin2_ckpt65x_20160612, ctrb_darwin2_ckpt66f_20170407, ctrb_darwin2_ckpt66a_20161020, ctrb_darwin2_ckpt66b_20161219, ctrb_darwin2_ckpt66j_20170815, ctrb_darwin2_ckpt65y_20160801, ctrb_darwin2_ckpt66c_20170121, ctrb_darwin2_ckpt65t_20160221, ctrb_darwin2_ckpt66o_20180209, ctrb_darwin2_ckpt66e_20170314, ctrb_darwin2_ckpt65u_20160315, ctrb_darwin2_ckpt65r_20151221, ctrb_darwin2_ckpt66i_20170718, ctrb_darwin2_ckpt66l_20171025, ctrb_darwin2_ckpt66m_20171213, HEAD
Changes since 1.3: +84 -24 lines
Ben Ward - some superficial structural changes allowing runs with no pfts
         - more significant structural and parameter changes to follow later

1 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F,v 1.3 2013/12/27 17:29:00 jahn 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 #ifdef ALLOW_QUOTA
11
12 c=============================================================
13 c subroutine quota_forcing
14 c step forward bio-chemical tracers in time
15 C==============================================================
16 SUBROUTINE QUOTA_FORCING(
17 U Ptr,
18 I bi,bj,imin,imax,jmin,jmax,
19 I myTime,myIter,myThid)
20 #include "SIZE.h"
21 #include "EEPARAMS.h"
22 #include "PARAMS.h"
23 #include "GRID.h"
24 #include "PTRACERS_SIZE.h"
25 #include "PTRACERS_PARAMS.h"
26 #include "GCHEM.h"
27 #include "QUOTA_SIZE.h"
28 #include "QUOTA.h"
29 #include "DARWIN_IO.h"
30 #include "DYNVARS.h"
31 #ifdef USE_QSW
32 #include "FFIELDS.h"
33 #endif
34
35 C === Global variables ===
36 c tracers
37 _RL Ptr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nDarwin)
38 INTEGER bi,bj,imin,imax,jmin,jmax
39 _RL myTime
40 INTEGER myIter
41 INTEGER myThid
42
43 C============== Local variables ============================================
44 c biomodel tracer arrays
45 _RL nutrient(iimax)
46 _RL biomass(iomax,npmax)
47 _RL orgmat(iomax-iChl,komax)
48 #ifdef FQUOTA
49 c iron partitioning
50 _RL freefe(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51 _RL freefu
52 _RL inputFel
53 #endif
54 c upstream arrays for sinking
55 _RL bioabove(iomax,npmax)
56 _RL orgabove(iomax-iChl,komax)
57 c some working variables
58 _RL sumpy
59 _RL sumpyup
60 c light variables
61 _RL PAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
62 _RL sfac(1-OLy:sNy+OLy)
63 _RL atten,lite
64 _RL newtime ! for sub-timestepping
65 _RL runtim ! time from tracer initialization
66 c
67 #ifdef ALLOW_DIAGNOSTICS
68 COJ for diagnostics
69 _RL PParr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
70 #endif
71 #ifdef ALLOW_TIMEAVE
72 #ifdef QUOTA_DIAG_LIMIT
73 _RL Rlim(iomax-iChl-1,npmax)
74 _RL Ilim(npmax)
75 _RL Tlim
76 _RL AP(iomax,npmax)
77 _RL HP(iomax,npmax)
78 #endif
79 #endif
80 c
81
82 c some local variables
83 _RL Tlocal
84 _RL Slocal
85 _RL PARlocal
86 _RL dzlocal
87 _RL dtplankton
88 _RL PP
89 c local tendencies
90 _RL dbiomass(iomax,npmax)
91 _RL dorgmat(iomax-iChl,komax)
92 _RL dnutrient(iimax)
93 _RL tmp
94
95 INTEGER bottom
96 INTEGER surface
97 INTEGER i,j,k,it,ktmp
98 INTEGER ii,io,jp,ko
99 INTEGER place
100 INTEGER debug
101 #ifdef ALLOW_DIAGNOSTICS
102 CHARACTER*8 diagname
103 #endif
104
105 c
106 c--------------------------------------------------
107 c initialise variables
108 DO j=1-OLy,sNy+OLy
109 DO i=1-OLx,sNx+OLx
110 do k=1,Nr
111 #ifdef FQUOTA
112 freefe(i,j,k) = 0.0 _d 0
113 # endif
114 PAR(i,j,k) = 0.0 _d 0
115 #ifdef ALLOW_DIAGNOSTICS
116 COJ for diagnostics
117 PParr(i,j,k) = 0. _d 0
118 #endif
119 enddo !k
120 ENDDO !i
121 ENDDO !j
122 c
123 c bio-chemical time loop
124 c--------------------------------------------------
125 DO it=1,nsubtime
126 c -------------------------------------------------
127 COJ cannot use dfloat because of adjoint
128 COJ division will be double precision anyway because of dTtracerLev
129 newtime=myTime-dTtracerLev(1)+
130 & float(it)*dTtracerLev(1)/float(nsubtime)
131 c print*,'it ',it,newtime,nsubtime,myTime
132 runtim=myTime-float(PTRACERS_Iter0)*dTtracerLev(1)
133
134 #ifdef FQUOTA
135 c determine iron partitioning - solve for free iron
136 call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
137 & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
138 & myIter, mythid)
139 #endif
140
141 c find light in each grid cell
142 c ---------------------------
143 c determine incident light
144 #ifndef READ_PAR
145 #ifdef USE_QSW
146 DO j=1-OLy,sNy+OLy
147 DO i=1-OLx,sNx+OLx
148 sur_par(i,j,bi,bj)=-parfrac*Qsw(i,j,bi,bj)*
149 & parconv*maskC(i,j,1,bi,bj)
150 ENDDO
151 ENDDO
152 #else
153 DO j=1-OLy,sNy+OLy
154 sfac(j)=0. _d 0
155 ENDDO
156 call darwin_insol(newTime,sfac,bj)
157 DO j=1-OLy,sNy+OLy
158 DO i=1-OLx,sNx+OLx
159 sur_par(i,j,bi,bj)=sfac(j)*maskC(i,j,1,bi,bj)/86400. _d 6
160 c if (i.eq.1.and.j.ge.1.and.j.le.sNy)
161 c & write(24,*) sur_par(i,j,bi,bj)
162 ENDDO
163 ENDDO
164 #endif
165 #endif
166
167 C.................................................................
168 C.................................................................
169
170
171 DO j=1,sNy
172 DO i=1,sNx
173 c surface PAR
174 c take ice coverage into account
175 #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
176 COJ ice coverage already taken into account by seaice package
177 lite=sur_par(i,j,bi,bj)
178 #else
179 #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
180 c if using Qsw and seaice, then ice fraction is already
181 c taken into account
182 lite=sur_par(i,j,bi,bj)
183 #else
184 lite=sur_par(i,j,bi,bj)*(1. _d 0-fice(i,j,bi,bj))
185 #endif
186 #endif
187 atten = 0. _d 0
188 sumpy = 0. _d 0
189 c
190 c FOR EACH LAYER ...
191 do k= 1, NR
192 if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
193 c ---------------------------------------------------------------------
194 c benw
195 c
196 c Fetch biomodel variables from ptr (ptracers)
197 c (making sure they are .ge. 0 - brute force)
198 c
199 c (set biomodel tendencies to zero, at the same time)
200 c
201 c *********************************************************************
202 place = 0
203 c Inorganic Nutrients
204 do ii=1,iimax
205 place = place + 1
206 c ambient nutrients for each element (1 to iimax)
207 nutrient(ii) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
208 dnutrient(ii) = 0. _d 0
209 enddo ! ii
210 c *********************************************************************
211 c Unicellular biomass (including chlorophyll biomass - for non-grazers)
212 do io=1,iomax
213 do jp=1,npmax
214 if (io.ne.iChlo.or.autotrophy(jp).gt.0. _d 0) then ! no grazer chlorophyll
215 place = place + 1
216 biomass(io,jp) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
217 ! biomasses above current layer for sinking
218 if (k.eq.1) then
219 bioabove(io,jp)=0. _d 0
220 endif
221 ! initialise biomass rate of change
222 dbiomass(io,jp) = 0. _d 0
223 else ! if grazer, fill chl biomass with zeros
224 biomass(io,jp) = 0. _d 0
225 endif
226 enddo ! jp
227 enddo
228 c *********************************************************************
229 c Organic matter
230 do io=1,iomax-iChl
231 do ko=1,komax
232 c mass of element x for all OM classes
233 place = place + 1
234 orgmat(io,ko) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
235 ! biomasses above current layer for sinking
236 if (k.eq.1) then
237 orgabove(io,ko) = 0. _d 0
238 endif
239 #ifdef SQUOTA
240 if (ko.and.1.and.io.eq.iSili) then
241 place = place - 1
242 orgmat(iSili,1) = 0. _d 0
243 orgabove(iSili,1) = 0. _d 0
244 endif
245 #endif
246 dorgmat(io,ko) = 0. _d 0
247 enddo ! ko
248 enddo ! io
249 c *********************************************************************
250 c
251 c ---------------------------------------------------------------------
252
253
254 c find local light for level k
255 sumpyup = sumpy
256 sumpy = 0. _d 0
257 do jp=1,npmax
258 #ifndef GEIDER
259 ! sum nitrogen biomass
260 sumpy = sumpy + biomass(iNitr,jp)
261 #else
262 ! sum chlorophyll
263 sumpy = sumpy + biomass(iChlo,jp)
264 #endif
265 enddo
266
267 atten= atten + (k_w + k_chl*sumpy)*5. _d -1*drF(k)
268 if (k.gt.1)then
269 atten = atten + (k_w+k_chl*sumpyup)*5. _d -1*drF(k-1)
270 endif
271 PAR(i,j,k) = lite*exp(-atten)
272 c
273 c Physical variables
274 PARlocal = PAR(i,j,k)
275 Tlocal = theta(i,j,k,bi,bj)
276 Slocal = salt(i,j,k,bi,bj)
277 c Free Iron
278 #ifdef FQUOTA
279 freefu = max(freefe(i,j,k),0. _d 0)
280 if (k.eq.1) then
281 inputFel = inputFe(i,j,bi,bj)
282 else
283 inputFel = 0. _d 0
284 endif
285 #endif
286 c Layer thickness
287 dzlocal = drF(k)*HFacC(i,j,k,bi,bj)
288 c
289 c set bottom=1.0 if the layer below is not ocean
290 ktmp=min(nR,k+1)
291 if(hFacC(i,j,ktmp,bi,bj).eq.0. _d 0.or.k.eq.Nr) then
292 bottom = 1
293 else
294 bottom = 0
295 endif
296 if (k.eq.1) then
297 surface = 1
298 else
299 surface = 0
300 endif
301
302 c set other arguments to zero
303 debug=0
304
305 if (debug.eq.7) print*,'Inorganic nutrients',nutrient
306 if (debug.eq.7) print*,'Plankton biomass', biomass
307 if (debug.eq.7) print*,'Organic nutrients',orgmat
308 if (debug.eq.8) print*,'k, PARlocal, dzlocal',
309 & k,PARlocal,dzlocal
310 c ---------------------------------------------------------------------
311 CALL QUOTA_PLANKTON(
312 I biomass, orgmat, nutrient,
313 O PP,
314 I bioabove,
315 I orgabove,
316 #ifdef FQUOTA
317 I freefu, inputFel,
318 #endif
319 #ifdef ALLOW_TIMEAVE
320 #ifdef QUOTA_DIAG_LIMIT
321 O AP, HP,
322 O Rlim, Ilim, Tlim,
323 #endif
324 #endif
325 I PARlocal, Tlocal, Slocal,
326 I bottom, surface, dzlocal,
327 O dbiomass, dorgmat, dnutrient,
328 I debug,
329 I runtim,
330 I MyThid)
331 c ---------------------------------------------------------------------
332 c
333 #ifdef RELAX_NUTS
334 if (darwin_relaxscale.gt.0. _d 0) then
335 !
336 IF ( darwin_NO3_relaxFile .NE. ' ' ) THEN
337 tmp=(Ptr(i,j,k,bi,bj,iNO3 )-no3_obs(i,j,k,bi,bj))
338 if (tmp.lt.0. _d 0) then
339 dnutrient(iNO3)=dnutrient(iNO3)
340 & -(tmp/darwin_relaxscale)
341 endif
342 ENDIF
343 #ifdef PQUOTA
344 IF ( darwin_PO4_relaxFile .NE. ' ' ) THEN
345 tmp=(Ptr(i,j,k,bi,bj,iPO4 )-po4_obs(i,j,k,bi,bj))
346 if (tmp.lt.0. _d 0) then
347 dnutrient(iPO4)=dnutrient(iPO4)
348 & -(tmp/darwin_relaxscale)
349 endif
350 ENDIF
351 #endif
352 #ifdef FQOUTA
353 IF ( darwin_Fet_relaxFile .NE. ' ' ) THEN
354 tmp=(Ptr(i,j,k,bi,bj,iFeT )-fet_obs(i,j,k,bi,bj))
355 if (tmp.lt.0. _d 0) then
356 dnutrient(iFeT)=dnutrient(iFeT)
357 & -(tmp/darwin_relaxscale)
358 endif
359 ENDIF
360 #endif
361 #ifdef SQUOTA
362 IF ( darwin_Si_relaxFile .NE. ' ' ) THEN
363 tmp=( Ptr(i,j,k,bi,bj,iSi )-si_obs(i,j,k,bi,bj))
364 if (tmp.lt.0. _d 0) then
365 dnutrient(iSi)=dnutrient(iSi)
366 & -(tmp/darwin_relaxscale)
367 endif
368 ENDIF
369 #endif
370 endif
371 #endif
372 c
373 #ifdef FQUOTA
374 #ifdef IRON_SED_SOURCE
375 c only above minimum depth (continental shelf)
376 if (rF(k).lt.depthfesed) then
377 c only if bottom layer
378 if (HFacC(i,j,k+1,bi,bj).eq.0. _d 0) then
379 #ifdef IRON_SED_SOURCE_VARIABLE
380 c calculate sink of POC into bottom layer
381 tmp=orgsink(2)*orgabove(iCarb,2)/dzlocal
382 c convert to dPOCl
383 dnutrient(iFeT) = dnutrient(iFeT)
384 & + fesedflux_pcm*tmp
385 #else
386 dnutrient(iFeT) = dnutrient(iFeT)
387 & + fesedflux/(drF(k)*hFacC(i,j,k,bi,bj))
388 #endif
389 endif
390 endif
391 #endif
392 #endif
393 c ---------------------------------------------------------------------
394 c save un-updated biomass as layer above
395 do io=1,iomax
396 do jp=1,npmax
397 bioabove(io,jp)=biomass(io,jp)
398 enddo
399 if (io.ne.iChlo) then
400 do ko=1,komax
401 orgabove(io,ko)=orgmat(io,ko)
402 enddo
403 endif
404 enddo
405 c ---------------------------------------------------------------------
406 c now update main tracer arrays
407 c for timestep dtplankton
408 dtplankton = dTtracerLev(k)/float(nsubtime)
409 cccccccccccccccccccccccccccccccccccccccccccccccccccc
410 place = 0
411 cccccccccccccccccccccccccccccccccccccccccccccccccccc
412 c Inorganic nutrients
413 do ii=1,iimax
414 place = place + 1
415 Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
416 & + dtplankton*dnutrient(ii)
417 enddo ! ii
418 cccccccccccccccccccccccccccccccccccccccccccccccccccc
419 c Biomass
420 do io=1,iomax
421 do jp=1,npmax
422 if (io.ne.iChlo.or.autotrophy(jp).gt.0. _d 0) then ! if not a grazer
423 place = place + 1
424 Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
425 & + dtplankton*dbiomass(io,jp)
426 endif
427 enddo ! jp
428 enddo ! io
429 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
430 c Organic matter
431 do io=1,iomax-iChl
432 do ko=1,komax
433 if (ko.ne.1.or.io.ne.iSili) then
434 place = place + 1
435 Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
436 & + dtplankton*dorgmat(io,ko)
437 endif
438 enddo ! ko
439 enddo ! io
440 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
441 c
442 #ifdef ALLOW_DIAGNOSTICS
443 COJ for diagnostics
444 PParr(i,j,k) = PP
445 #endif /* ALLOW_DIAGNOSTICS */
446
447 #ifdef ALLOW_TIMEAVE
448 PPave(i,j,k,bi,bj) = PPave(i,j,k,bi,bj)
449 & + PP * dtplankton
450 PARave(i,j,k,bi,bj) = PARave(i,j,k,bi,bj)
451 & + PARlocal * dtplankton
452 c
453 #ifdef QUOTA_DIAG_LIMIT
454 do jp=1,npmax
455 ! carbon
456 AP_C_ave(i,j,k,bi,bj,jp) = AP_C_ave(i,j,k,bi,bj,jp)
457 & + AP(iCarb,jp) * dtplankton
458 HP_C_ave(i,j,k,bi,bj,jp) = HP_C_ave(i,j,k,bi,bj,jp)
459 & + HP(iCarb,jp) * dtplankton
460 ! nitrogen
461 AP_N_ave(i,j,k,bi,bj,jp) = AP_N_ave(i,j,k,bi,bj,jp)
462 & + AP(iNitr,jp) * dtplankton
463 HP_N_ave(i,j,k,bi,bj,jp) = HP_N_ave(i,j,k,bi,bj,jp)
464 & + HP(iNitr,jp) * dtplankton
465 Nlimave(i,j,k,bi,bj,jp) = Nlimave(i,j,k,bi,bj,jp)
466 & + Rlim(iNitr-1,jp) * dtplankton
467 ! phosphorus
468 #ifdef PQUOTA
469 AP_P_ave(i,j,k,bi,bj,jp) = AP_P_ave(i,j,k,bi,bj,jp)
470 & + AP(iPhos,jp) * dtplankton
471 HP_P_ave(i,j,k,bi,bj,jp) = HP_P_ave(i,j,k,bi,bj,jp)
472 & + HP(iPhos,jp) * dtplankton
473 Plimave(i,j,k,bi,bj,jp) = Plimave(i,j,k,bi,bj,jp)
474 & + Rlim(iPhos-1,jp) * dtplankton
475 #endif
476 ! iron
477 #ifdef FQUOTA
478 AP_F_ave(i,j,k,bi,bj,jp) = AP_F_ave(i,j,k,bi,bj,jp)
479 & + AP(iIron,jp) * dtplankton
480 HP_F_ave(i,j,k,bi,bj,jp) = HP_F_ave(i,j,k,bi,bj,jp)
481 & + HP(iIron,jp) * dtplankton
482 Flimave(i,j,k,bi,bj,jp) = Flimave(i,j,k,bi,bj,jp)
483 & + Rlim(iIron-1,jp) * dtplankton
484 #endif
485 ! light
486 Ilimave(i,j,k,bi,bj,jp) = Ilimave(i,j,k,bi,bj,jp)
487 & + Ilim(jp) * dtplankton
488 enddo
489 Tlimave(i,j,k,bi,bj) = Tlimave(i,j,k,bi,bj)
490 & + Tlim * dtplankton
491 #endif
492 #endif
493 endif
494 c end if hFac>0
495 enddo ! k
496 c end layer loop
497 c
498 ENDDO ! i
499 ENDDO ! j
500 c
501 c
502 COJ fill diagnostics
503 #ifdef ALLOW_DIAGNOSTICS
504 IF ( useDiagnostics ) THEN
505 diagname = 'PP '
506 CALL DIAGNOSTICS_FILL( PParr(1-Olx,1-Oly,1), diagname,
507 & 0,Nr,2,bi,bj,myThid )
508 ENDIF
509 #endif
510 COJ
511
512 #ifdef FQUOTA
513 c determine iron partitioning - solve for free iron
514 call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
515 & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
516 & myIter, mythid)
517 #endif
518
519 c
520 #ifdef ALLOW_TIMEAVE
521 c save averages
522 dar_timeave(bi,bj) = dar_timeave(bi,bj) + dtplankton
523 #endif
524 c
525 c -----------------------------------------------------
526 ENDDO ! it
527 c -----------------------------------------------------
528 c end of bio-chemical time loop
529 c
530 RETURN
531 END
532
533 #endif /*ALLOW_QUOTA*/
534 #endif /*ALLOW_DARWIN*/
535 #endif /*ALLOW_PTRACERS*/
536
537 C============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22