/[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.4 - (hide 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 benw 1.4 C $Header: /u/gcmpack/MITgcm_contrib/darwin2/pkg/quota/quota_forcing.F,v 1.3 2013/12/27 17:29:00 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 benw 1.4 c upstream arrays for sinking
55 jahn 1.1 _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 benw 1.2 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 benw 1.4 _RL Rlim(iomax-iChl-1,npmax)
74 benw 1.2 _RL Ilim(npmax)
75     _RL Tlim
76 benw 1.4 _RL AP(iomax,npmax)
77     _RL HP(iomax,npmax)
78 benw 1.2 #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.4 INTEGER i,j,k,it,ktmp
98     INTEGER ii,io,jp,ko
99 jahn 1.1 INTEGER place
100     INTEGER debug
101 benw 1.4 #ifdef ALLOW_DIAGNOSTICS
102 jahn 1.1 CHARACTER*8 diagname
103 benw 1.4 #endif
104 jahn 1.1
105     c
106     c--------------------------------------------------
107 benw 1.2 c initialise variables
108 jahn 1.1 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 benw 1.2 #ifdef ALLOW_DIAGNOSTICS
116     COJ for diagnostics
117     PParr(i,j,k) = 0. _d 0
118 jahn 1.1 #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 benw 1.2 c Unicellular biomass (including chlorophyll biomass - for non-grazers)
212 jahn 1.1 do io=1,iomax
213     do jp=1,npmax
214 benw 1.4 if (io.ne.iChlo.or.autotrophy(jp).gt.0. _d 0) then ! no grazer chlorophyll
215 benw 1.2 place = place + 1
216     biomass(io,jp) = max(Ptr(i,j,k,bi,bj,place),0. _d 0)
217 jahn 1.1 ! biomasses above current layer for sinking
218 benw 1.2 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 jahn 1.1 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 benw 1.4 I bioabove,
315 jahn 1.1 I orgabove,
316     #ifdef FQUOTA
317     I freefu, inputFel,
318     #endif
319 benw 1.2 #ifdef ALLOW_TIMEAVE
320     #ifdef QUOTA_DIAG_LIMIT
321 benw 1.4 O AP, HP,
322     O Rlim, Ilim, Tlim,
323 benw 1.2 #endif
324     #endif
325 jahn 1.1 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 benw 1.4 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 benw 1.2 #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 jahn 1.1 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 benw 1.4 if (io.ne.iChlo.or.autotrophy(jp).gt.0. _d 0) then ! if not a grazer
423 benw 1.2 place = place + 1
424     Ptr(i,j,k,bi,bj,place) = Ptr(i,j,k,bi,bj,place)
425     & + dtplankton*dbiomass(io,jp)
426 jahn 1.1 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 benw 1.2 #ifdef ALLOW_DIAGNOSTICS
443     COJ for diagnostics
444     PParr(i,j,k) = PP
445     #endif /* ALLOW_DIAGNOSTICS */
446    
447 jahn 1.1 #ifdef ALLOW_TIMEAVE
448 benw 1.2 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 jahn 1.1 c
453 benw 1.2 #ifdef QUOTA_DIAG_LIMIT
454     do jp=1,npmax
455 benw 1.4 ! 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 benw 1.2 Nlimave(i,j,k,bi,bj,jp) = Nlimave(i,j,k,bi,bj,jp)
466 benw 1.4 & + 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 benw 1.2 Flimave(i,j,k,bi,bj,jp) = Flimave(i,j,k,bi,bj,jp)
483 benw 1.4 & + Rlim(iIron-1,jp) * dtplankton
484     #endif
485     ! light
486 benw 1.2 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 jahn 1.1 #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 benw 1.2 diagname = 'PP '
506     CALL DIAGNOSTICS_FILL( PParr(1-Olx,1-Oly,1), diagname,
507 jahn 1.1 & 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 jahn 1.3 dar_timeave(bi,bj) = dar_timeave(bi,bj) + dtplankton
523 jahn 1.1 #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