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

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

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


Revision 1.2 - (hide annotations) (download)
Mon Jul 2 09:47:43 2012 UTC (13 years, 1 month ago) by benw
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64k_20130723, ctrb_darwin2_ckpt64h_20130528, ctrb_darwin2_ckpt64m_20130820, ctrb_darwin2_ckpt64r_20131210, ctrb_darwin2_ckpt64f_20130405, ctrb_darwin2_ckpt64a_20121116, ctrb_darwin2_ckpt64n_20130826, ctrb_darwin2_ckpt64o_20131024, ctrb_darwin2_ckpt64i_20130622, ctrb_darwin2_ckpt64e_20130305, ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt63r_20120817, ctrb_darwin2_ckpt64g_20130503, ctrb_darwin2_ckpt64l_20130806, ctrb_darwin2_ckpt64c_20130120, ctrb_darwin2_ckpt64j_20130704, ctrb_darwin2_ckpt63p_20120707, ctrb_darwin2_ckpt64p_20131118, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt64b_20121224, ctrb_darwin2_ckpt64d_20130219, ctrb_darwin2_ckpt64_20121012, ctrb_darwin2_ckpt64q_20131118, ctrb_darwin2_ckpt64p_20131024
Changes since 1.1: +91 -76 lines
Removed obsolete diversity diagnostics
Added N,Fe,T,I limitation factors
Added option to not carry zero chlorophyll tracer for zooplankton (new option iPhoto in QUOTA_SIZE.h)

1 benw 1.2 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F,v 1.1 2011/04/13 18:56:26 jahn Exp $
2 jahn 1.1 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/swimming
55     _RL bioabove(iomax,npmax)
56     _RL biobelow(iomax,npmax)
57     _RL orgabove(iomax-iChl,komax)
58     c some working variables
59     _RL sumpy
60     _RL sumpyup
61     c light variables
62     _RL PAR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
63     _RL sfac(1-OLy:sNy+OLy)
64     _RL atten,lite
65     _RL newtime ! for sub-timestepping
66     _RL runtim ! time from tracer initialization
67 benw 1.2 c
68     #ifdef ALLOW_DIAGNOSTICS
69     COJ for diagnostics
70     _RL PParr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
71     #endif
72     #ifdef ALLOW_TIMEAVE
73     #ifdef QUOTA_DIAG_LIMIT
74     _RL Nlim(npmax)
75     _RL Flim(npmax)
76     _RL Ilim(npmax)
77     _RL Tlim
78     #endif
79 jahn 1.1 #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 benw 1.2 INTEGER i,j,k,it,itmp,ktmp
98 jahn 1.1 INTEGER ii,io,jp,ko, jp2, jpsave
99     INTEGER place
100     INTEGER debug
101     CHARACTER*8 diagname
102    
103     c
104     c--------------------------------------------------
105 benw 1.2 c initialise variables
106 jahn 1.1 DO j=1-OLy,sNy+OLy
107     DO i=1-OLx,sNx+OLx
108     do k=1,Nr
109     #ifdef FQUOTA
110     freefe(i,j,k) = 0.0 _d 0
111     # endif
112     PAR(i,j,k) = 0.0 _d 0
113 benw 1.2 #ifdef ALLOW_DIAGNOSTICS
114     COJ for diagnostics
115     PParr(i,j,k) = 0. _d 0
116 jahn 1.1 #endif
117     enddo !k
118     ENDDO !i
119     ENDDO !j
120     c
121     c bio-chemical time loop
122     c--------------------------------------------------
123     DO it=1,nsubtime
124     c -------------------------------------------------
125     COJ cannot use dfloat because of adjoint
126     COJ division will be double precision anyway because of dTtracerLev
127     newtime=myTime-dTtracerLev(1)+
128     & float(it)*dTtracerLev(1)/float(nsubtime)
129     c print*,'it ',it,newtime,nsubtime,myTime
130     runtim=myTime-float(PTRACERS_Iter0)*dTtracerLev(1)
131    
132     #ifdef FQUOTA
133     c determine iron partitioning - solve for free iron
134     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
135     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
136     & myIter, mythid)
137     #endif
138    
139     c find light in each grid cell
140     c ---------------------------
141     c determine incident light
142     #ifndef READ_PAR
143     #ifdef USE_QSW
144     DO j=1-OLy,sNy+OLy
145     DO i=1-OLx,sNx+OLx
146     sur_par(i,j,bi,bj)=-parfrac*Qsw(i,j,bi,bj)*
147     & parconv*maskC(i,j,1,bi,bj)
148     ENDDO
149     ENDDO
150     #else
151     DO j=1-OLy,sNy+OLy
152     sfac(j)=0. _d 0
153     ENDDO
154     call darwin_insol(newTime,sfac,bj)
155     DO j=1-OLy,sNy+OLy
156     DO i=1-OLx,sNx+OLx
157     sur_par(i,j,bi,bj)=sfac(j)*maskC(i,j,1,bi,bj)/86400. _d 6
158     c if (i.eq.1.and.j.ge.1.and.j.le.sNy)
159     c & write(24,*) sur_par(i,j,bi,bj)
160     ENDDO
161     ENDDO
162     #endif
163     #endif
164    
165     C.................................................................
166     C.................................................................
167    
168    
169     DO j=1,sNy
170     DO i=1,sNx
171     c surface PAR
172     c take ice coverage into account
173     #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
174     COJ ice coverage already taken into account by seaice package
175     lite=sur_par(i,j,bi,bj)
176     #else
177     #if (defined (ALLOW_SEAICE) && defined (USE_QSW))
178     c if using Qsw and seaice, then ice fraction is already
179     c taken into account
180     lite=sur_par(i,j,bi,bj)
181     #else
182     lite=sur_par(i,j,bi,bj)*(1. _d 0-fice(i,j,bi,bj))
183     #endif
184     #endif
185     atten = 0. _d 0
186     sumpy = 0. _d 0
187     c
188     c FOR EACH LAYER ...
189     do k= 1, NR
190     if (HFacC(i,j,k,bi,bj).gt.0. _d 0) then
191     c ---------------------------------------------------------------------
192     c benw
193     c
194     c Fetch biomodel variables from ptr (ptracers)
195     c (making sure they are .ge. 0 - brute force)
196     c
197     c (set biomodel tendencies to zero, at the same time)
198     c
199     c *********************************************************************
200     place = 0
201     c Inorganic Nutrients
202     do ii=1,iimax
203     place = place + 1
204     c ambient nutrients for each element (1 to iimax)
205     nutrient(ii) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
206     dnutrient(ii) = 0. _d 0
207     enddo ! ii
208     c *********************************************************************
209 benw 1.2 c Unicellular biomass (including chlorophyll biomass - for non-grazers)
210 jahn 1.1 do io=1,iomax
211     do jp=1,npmax
212 benw 1.2 if (io.ne.iChlo.or.pft(jp).ne.6) then ! no grazer chlorophyll
213     place = place + 1
214     biomass(io,jp) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
215 jahn 1.1 ! biomasses above current layer for sinking
216 benw 1.2 if (k.eq.1) then
217     bioabove(io,jp)=0. _d 0
218     endif
219 jahn 1.1 ! biomasses below current layer for swimming
220 benw 1.2 if (k.eq.Nr) then
221     biobelow(io,jp)=0. _d 0
222     elseif (hFacC(i,j,k+1,bi,bj).eq.0. _d 0) then
223     biobelow(io,jp)=0. _d 0
224     else
225     biobelow(io,jp)=max(Ptr(i,j,k+1,bi,bj,place),0. _d 0)
226     endif
227     ! initialise biomass rate of change
228     dbiomass(io,jp) = 0. _d 0
229     else ! if grazer, fill chl biomass with zeros
230     biomass(io,jp) = 0. _d 0
231 jahn 1.1 endif
232     enddo ! jp
233     enddo
234     c *********************************************************************
235     c Organic matter
236     do io=1,iomax-iChl
237     do ko=1,komax
238     c mass of element x for all OM classes
239     place = place + 1
240     orgmat(io,ko) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
241     ! biomasses above current layer for sinking
242     if (k.eq.1) then
243     orgabove(io,ko) = 0. _d 0
244     endif
245     #ifdef SQUOTA
246     if (ko.and.1.and.io.eq.iSili) then
247     place = place - 1
248     orgmat(iSili,1) = 0. _d 0
249     orgabove(iSili,1) = 0. _d 0
250     endif
251     #endif
252     dorgmat(io,ko) = 0. _d 0
253     enddo ! ko
254     enddo ! io
255     c *********************************************************************
256     c
257     c ---------------------------------------------------------------------
258    
259    
260     c find local light for level k
261     sumpyup = sumpy
262     sumpy = 0. _d 0
263     do jp=1,npmax
264     #ifndef GEIDER
265     ! sum nitrogen biomass
266     sumpy = sumpy + biomass(iNitr,jp)
267     #else
268     ! sum chlorophyll
269     sumpy = sumpy + biomass(iChlo,jp)
270     #endif
271     enddo
272    
273     atten= atten + (k_w + k_chl*sumpy)*5. _d -1*drF(k)
274     if (k.gt.1)then
275     atten = atten + (k_w+k_chl*sumpyup)*5. _d -1*drF(k-1)
276     endif
277     PAR(i,j,k) = lite*exp(-atten)
278     c
279     c Physical variables
280     PARlocal = PAR(i,j,k)
281     Tlocal = theta(i,j,k,bi,bj)
282     Slocal = salt(i,j,k,bi,bj)
283     c Free Iron
284     #ifdef FQUOTA
285     freefu = max(freefe(i,j,k),0. _d 0)
286     if (k.eq.1) then
287     inputFel = inputFe(i,j,bi,bj)
288     else
289     inputFel = 0. _d 0
290     endif
291     #endif
292     c Layer thickness
293     dzlocal = drF(k)*HFacC(i,j,k,bi,bj)
294     c
295     c set bottom=1.0 if the layer below is not ocean
296     ktmp=min(nR,k+1)
297     if(hFacC(i,j,ktmp,bi,bj).eq.0. _d 0.or.k.eq.Nr) then
298     bottom = 1
299     else
300     bottom = 0
301     endif
302     if (k.eq.1) then
303     surface = 1
304     else
305     surface = 0
306     endif
307    
308     c set other arguments to zero
309     debug=0
310    
311     if (debug.eq.7) print*,'Inorganic nutrients',nutrient
312     if (debug.eq.7) print*,'Plankton biomass', biomass
313     if (debug.eq.7) print*,'Organic nutrients',orgmat
314     if (debug.eq.8) print*,'k, PARlocal, dzlocal',
315     & k,PARlocal,dzlocal
316     c ---------------------------------------------------------------------
317     CALL QUOTA_PLANKTON(
318     I biomass, orgmat, nutrient,
319     O PP,
320     I bioabove, biobelow,
321     I orgabove,
322     #ifdef FQUOTA
323     I freefu, inputFel,
324     #endif
325 benw 1.2 #ifdef ALLOW_TIMEAVE
326     #ifdef QUOTA_DIAG_LIMIT
327     O Nlim, Flim, Ilim, Tlim,
328     #endif
329     #endif
330 jahn 1.1 I PARlocal, Tlocal, Slocal,
331     I bottom, surface, dzlocal,
332     O dbiomass, dorgmat, dnutrient,
333     I debug,
334     I runtim,
335     I MyThid)
336     c ---------------------------------------------------------------------
337 benw 1.2 #ifdef FQUOTA
338     #ifdef IRON_SED_SOURCE
339     c only above minimum depth (continental shelf)
340     if (rF(k).lt.depthfesed) then
341     c only if bottom layer
342     if (HFacC(i,j,k+1,bi,bj).eq.0. _d 0) then
343     #ifdef IRON_SED_SOURCE_VARIABLE
344     c calculate sink of POC into bottom layer
345     tmp=orgsink(2)*orgabove(iCarb,2)/dzlocal
346     c convert to dPOCl
347     dnutrient(iFeT) = dnutrient(iFeT)
348     & + fesedflux_pcm*tmp
349     #else
350     dnutrient(iFeT) = dnutrient(iFeT)
351     & + fesedflux/(drF(k)*hFacC(i,j,k,bi,bj))
352     #endif
353     endif
354     endif
355     #endif
356     #endif
357     c ---------------------------------------------------------------------
358 jahn 1.1 c save un-updated biomass as layer above
359     do io=1,iomax
360     do jp=1,npmax
361     bioabove(io,jp)=biomass(io,jp)
362     enddo
363     if (io.ne.iChlo) then
364     do ko=1,komax
365     orgabove(io,ko)=orgmat(io,ko)
366     enddo
367     endif
368     enddo
369     c ---------------------------------------------------------------------
370     c now update main tracer arrays
371     c for timestep dtplankton
372     dtplankton = dTtracerLev(k)/float(nsubtime)
373     cccccccccccccccccccccccccccccccccccccccccccccccccccc
374     place = 0
375     cccccccccccccccccccccccccccccccccccccccccccccccccccc
376     c Inorganic nutrients
377     do ii=1,iimax
378     place = place + 1
379     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
380     & + dtplankton*dnutrient(ii)
381     enddo ! ii
382     cccccccccccccccccccccccccccccccccccccccccccccccccccc
383     c Biomass
384     do io=1,iomax
385     do jp=1,npmax
386 benw 1.2 if (io.ne.iChlo.or.pft(jp).ne.6) then ! if not a grazer
387     place = place + 1
388     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
389     & + dtplankton*dbiomass(io,jp)
390     if (pft(jp).eq.6.and.io.eq.iChlo) then
391     Ptr(i,j,k,bi,bj,place) = 0. _d 0
392     endif
393 jahn 1.1 endif
394     enddo ! jp
395     enddo ! io
396     ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
397     c Organic matter
398     do io=1,iomax-iChl
399     do ko=1,komax
400     if (ko.ne.1.or.io.ne.iSili) then
401     place = place + 1
402     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
403     & + dtplankton*dorgmat(io,ko)
404     endif
405     enddo ! ko
406     enddo ! io
407     ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
408     c
409 benw 1.2 #ifdef ALLOW_DIAGNOSTICS
410     COJ for diagnostics
411     PParr(i,j,k) = PP
412     #endif /* ALLOW_DIAGNOSTICS */
413    
414 jahn 1.1 #ifdef ALLOW_TIMEAVE
415 benw 1.2 PPave(i,j,k,bi,bj) = PPave(i,j,k,bi,bj)
416     & + PP * dtplankton
417     PARave(i,j,k,bi,bj) = PARave(i,j,k,bi,bj)
418     & + PARlocal * dtplankton
419 jahn 1.1 c
420 benw 1.2 #ifdef QUOTA_DIAG_LIMIT
421     do jp=1,npmax
422     Nlimave(i,j,k,bi,bj,jp) = Nlimave(i,j,k,bi,bj,jp)
423     & + Nlim(jp) * dtplankton
424     Flimave(i,j,k,bi,bj,jp) = Flimave(i,j,k,bi,bj,jp)
425     & + Flim(jp) * dtplankton
426     Ilimave(i,j,k,bi,bj,jp) = Ilimave(i,j,k,bi,bj,jp)
427     & + Ilim(jp) * dtplankton
428     enddo
429     Tlimave(i,j,k,bi,bj) = Tlimave(i,j,k,bi,bj)
430     & + Tlim * dtplankton
431 jahn 1.1 #endif
432     #endif
433     endif
434     c end if hFac>0
435     enddo ! k
436     c end layer loop
437     c
438     ENDDO ! i
439     ENDDO ! j
440     c
441     c
442     COJ fill diagnostics
443     #ifdef ALLOW_DIAGNOSTICS
444     IF ( useDiagnostics ) THEN
445 benw 1.2 diagname = 'PP '
446     CALL DIAGNOSTICS_FILL( PParr(1-Olx,1-Oly,1), diagname,
447 jahn 1.1 & 0,Nr,2,bi,bj,myThid )
448     ENDIF
449     #endif
450     COJ
451    
452     #ifdef FQUOTA
453     c determine iron partitioning - solve for free iron
454     call darwin_fe_chem(bi,bj,iMin,iMax,jMin,jMax,
455     & Ptr(1-OLx,1-OLy,1,bi,bj,iFeT), freefe,
456     & myIter, mythid)
457     #endif
458    
459     c
460     #ifdef ALLOW_TIMEAVE
461     c save averages
462     do k=1,nR
463     dar_timeave(bi,bj,k) = dar_timeave(bi,bj,k)
464     & + dtplankton
465     enddo
466     #endif
467     c
468     c -----------------------------------------------------
469     ENDDO ! it
470     c -----------------------------------------------------
471     c end of bio-chemical time loop
472     c
473     RETURN
474     END
475    
476     #endif /*ALLOW_QUOTA*/
477     #endif /*ALLOW_DARWIN*/
478     #endif /*ALLOW_PTRACERS*/
479    
480     C============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22