New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
phytoplankton.F90 in branches/NERC/dev_r5518_GO6_MEDUSA_conserv/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/NERC/dev_r5518_GO6_MEDUSA_conserv/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90 @ 8489

Last change on this file since 8489 was 8489, checked in by jpalmier, 7 years ago

JPALM -- gmed ticket #346 : improve MEDUSA conservation -- import BBL bug fix from NEMO ticket #1932 to GO6 branch

File size: 21.9 KB
Line 
1MODULE phytoplankton_mod
2   !!======================================================================
3   !!                         ***  MODULE phytoplankton_mod  ***
4   !! Calculates the phytoplankton growth
5   !!======================================================================
6   !! History :
7   !!   -   ! 2017-04 (M. Stringer)        Code taken from trcbio_medusa.F90
8   !!   -   ! 2017-08 (A. Yool)            Mean mixed layer chlorophyll
9   !!----------------------------------------------------------------------
10#if defined key_medusa
11   !!----------------------------------------------------------------------
12   !!                                                   MEDUSA bio-model
13   !!----------------------------------------------------------------------
14
15   IMPLICIT NONE
16   PRIVATE
17     
18   PUBLIC   phytoplankton        ! Called in plankton.F90
19
20   !!----------------------------------------------------------------------
21   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
22   !! $Id$
23   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
24   !!----------------------------------------------------------------------
25
26CONTAINS
27
28   SUBROUTINE phytoplankton( jk )
29      !!---------------------------------------------------------------------
30      !!                     ***  ROUTINE phytoplankton  ***
31      !! This called from PLANKTON and calculates the phytoplankton
32      !! growth.
33      !!----------------------------------------------------------------------
34      USE bio_medusa_mod,    ONLY: fdep1, ffld, ffln2,                   &
35                                   fjlim_pd, fjlim_pn,                   &
36                                   fnld, fnln,                           &
37                                   fprd, fprd_ml, fprds,                 &
38                                   fprn, fprn_ml, frd, frn,              &
39                                   fsin, fsld, fsld2, fthetad, fthetan,  & 
40                                   ftot_det, ftot_dtc, ftot_pd,          &
41                                   ftot_pn, ftot_zme, ftot_zmi,          &
42                                   fun_Q10, fun_T, idf, idfval,          &
43                                   zchd, zchn, zdet, zdin, zdtc,         &
44                                   zfer, zpds, zphd, zphn, zsil,         &
45                                   zzme, zzmi, fchl_ml
46      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
47      USE in_out_manager,    ONLY: lwp, numout
48      USE oce,               ONLY: tsn
49      USE par_kind,          ONLY: wp
50      USE par_oce,           ONLY: jp_tem, jpi, jpim1, jpj, jpjm1
51      USE phycst,            ONLY: rsmall
52      USE sms_medusa,        ONLY: jliebig, jphy, jq10,                  &
53                                   xald, xaln, xfld, xfln,               &
54                                   xnld, xnln, xnsi0, xpar,              &
55                                   xsin0, xsld, xthetam, xthetamd, xuif, &
56                                   xvpd, xvpn, xxi
57      USE zdfmxl,            ONLY: hmld
58
59   !!* Substitution
60#  include "domzgr_substitute.h90"
61
62      !! Level
63      INTEGER, INTENT( in ) :: jk
64
65      INTEGER :: ji, jj
66
67      REAL(wp), DIMENSION(jpi,jpj) :: faln, fchn, fjln
68      REAL(wp), DIMENSION(jpi,jpj) :: fald, fchd, fjld
69      REAL(wp)                     :: fchn1, fchd1
70      !! AXY (03/02/11): add in Liebig terms
71      REAL(wp)                     :: fpnlim, fpdlim
72      !! AXY (16/07/09): add in Eppley curve functionality
73      REAL(wp)                     :: xvpnT,xvpdT
74      !! silicon cycle
75      REAL(wp)                     :: fnsi
76
77      REAL(wp)                     :: fsin1, fnsi1, fnsi2
78      REAL(wp)                     :: fq0
79
80      DO jj = 2,jpjm1
81         DO ji = 2,jpim1
82            !! OPEN wet point IF..THEN loop
83            if (tmask(ji,jj,jk) == 1) then
84               !!----------------------------------------------------------
85               !! Chlorophyll calculations
86               !!----------------------------------------------------------
87               !!
88               !! non-diatoms
89          if (zphn(ji,jj).GT.rsmall) then
90                  fthetan(ji,jj) = max(tiny(zchn(ji,jj)),                    &
91                                       (zchn(ji,jj) * xxi) /                 &
92                                       (zphn(ji,jj) + tiny(zphn(ji,jj))))
93                  faln(ji,jj)    = xaln * fthetan(ji,jj)
94               else
95                  fthetan(ji,jj) = 0.
96                  faln(ji,jj)    = 0.
97               endif
98               !!
99               !! diatoms
100          if (zphd(ji,jj).GT.rsmall) then
101                  fthetad(ji,jj) = max(tiny(zchd(ji,jj)),                   &
102                                       (zchd(ji,jj) * xxi) /                &
103                                       (zphd(ji,jj) + tiny(zphd(ji,jj))))
104                  fald(ji,jj)    = xald * fthetad(ji,jj)
105               else
106                  fthetad(ji,jj) = 0.
107                  fald(ji,jj)    = 0.
108               endif
109
110# if defined key_debug_medusa
111               !! report biological calculations
112               if (idf.eq.1.AND.idfval.eq.1) then
113                  IF (lwp) write (numout,*) '------------------------------'
114                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln(ji,jj)
115                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald(ji,jj)
116               endif
117# endif
118            ENDIF
119         ENDDO
120      ENDDO
121
122      DO jj = 2,jpjm1
123         DO ji = 2,jpim1
124            if (tmask(ji,jj,jk) == 1) then
125               !!----------------------------------------------------------
126               !! Phytoplankton light limitation
127               !!----------------------------------------------------------
128               !!
129               !! It is assumed xpar is the depth-averaged (vertical layer) PAR
130               !! Light limitation (check self-shading) in W/m2
131               !!
132               !! Note that there is no temperature dependence in phytoplankton
133               !! growth rate or any other function.
134               !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to
135               !! avoid NaNs in case of Phy==0. 
136               !!
137               !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and
138               !! non-diat:
139               !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012
140               !!
141               !! AXY (16/07/09)
142               !! temperature for new Eppley style phytoplankton growth
143               fun_T(ji,jj)   = 1.066**(1.0 * tsn(ji,jj,jk,jp_tem))
144               !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
145               !!                 phytoplankton growth; remin. unaffected
146               fun_Q10(ji,jj) = jq10**((tsn(ji,jj,jk,jp_tem) - 0.0) / 10.0)
147               if (jphy.eq.1) then
148                  xvpnT = xvpn * fun_T(ji,jj)
149                  xvpdT = xvpd * fun_T(ji,jj)
150               elseif (jphy.eq.2) then
151                  xvpnT = xvpn * fun_Q10(ji,jj)
152                  xvpdT = xvpd * fun_Q10(ji,jj)
153               else
154                  xvpnT = xvpn
155                  xvpdT = xvpd
156               endif
157               !!
158               !! non-diatoms
159               fchn1 = (xvpnT * xvpnT) +                                     &
160                       (faln(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk) *         &
161                        xpar(ji,jj,jk))
162               if (fchn1.GT.rsmall) then
163                  fchn(ji,jj) = xvpnT / (sqrt(fchn1) + tiny(fchn1))
164               else
165                  fchn(ji,jj) = 0.
166               endif
167               !! non-diatom J term
168               fjln(ji,jj)     = fchn(ji,jj) * faln(ji,jj) * xpar(ji,jj,jk)
169               fjlim_pn(ji,jj) = fjln(ji,jj) / xvpnT
170               !!
171               !! diatoms
172               fchd1 = (xvpdT * xvpdT) +                                     &
173                       (fald(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk) *         &
174                        xpar(ji,jj,jk))
175               if (fchd1.GT.rsmall) then
176                  fchd(ji,jj) = xvpdT / (sqrt(fchd1) + tiny(fchd1))
177               else
178                  fchd(ji,jj) = 0.
179               endif
180               !! diatom J term
181               fjld(ji,jj)    = fchd(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk)
182               fjlim_pd(ji,jj) = fjld(ji,jj) / xvpdT
183     
184# if defined key_debug_medusa
185               !! report phytoplankton light limitation
186               if (idf.eq.1.AND.idfval.eq.1) then
187                  IF (lwp) write (numout,*) '------------------------------'
188                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn(ji,jj)
189                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd(ji,jj)
190                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln(ji,jj)
191                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld(ji,jj)
192               endif
193# endif
194            ENDIF
195         ENDDO
196      ENDDO
197
198      DO jj = 2,jpjm1
199         DO ji = 2,jpim1
200            if (tmask(ji,jj,jk) == 1) then
201               !!----------------------------------------------------------
202               !! Phytoplankton nutrient limitation
203               !!----------------------------------------------------------
204               !!
205               !! non-diatoms (N, Fe).
206               !! non-diatom Qn term
207               fnln(ji,jj)  = zdin(ji,jj) / (zdin(ji,jj) + xnln)
208               !! non-diatom Qf term
209               ffln2(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfln)
210               !!
211               !! diatoms (N, Si, Fe).
212               !! diatom Qn term
213               fnld(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnld)
214               !! diatom Qs term
215               fsld(ji,jj) = zsil(ji,jj) / (zsil(ji,jj) + xsld)
216               !! diatom Qf term
217               ffld(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfld)
218
219# if defined key_debug_medusa
220               !! report phytoplankton nutrient limitation
221               if (idf.eq.1.AND.idfval.eq.1) then
222                  IF (lwp) write (numout,*) '------------------------------'
223                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln(ji,jj)
224                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld(ji,jj)
225                  IF (lwp) write (numout,*) 'ffln2(',jk,') = ', ffln2(ji,jj)
226                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld(ji,jj)
227                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld(ji,jj)
228               endif
229# endif
230            ENDIF
231         ENDDO
232      ENDDO
233
234      DO jj = 2,jpjm1
235         DO ji = 2,jpim1
236            if (tmask(ji,jj,jk) == 1) then
237               !!----------------------------------------------------------
238               !! Primary production (non-diatoms)
239               !! (note: still needs multiplying by phytoplankton
240               !! concentration)
241               !!----------------------------------------------------------
242               !!
243               if (jliebig .eq. 0) then
244                  !! multiplicative nutrient limitation
245                  fpnlim = fnln(ji,jj) * ffln2(ji,jj)
246               elseif (jliebig .eq. 1) then
247                  !! Liebig Law (= most limiting) nutrient limitation
248                  fpnlim = min(fnln(ji,jj), ffln2(ji,jj))
249               endif
250               fprn(ji,jj) = fjln(ji,jj) * fpnlim
251            ENDIF
252         ENDDO
253      ENDDO
254
255      DO jj = 2,jpjm1
256         DO ji = 2,jpim1
257            if (tmask(ji,jj,jk) == 1) then
258               !!----------------------------------------------------------
259               !! Primary production (diatoms)
260               !! (note: still needs multiplying by phytoplankton
261               !! concentration)
262               !!
263               !! Production here is split between nitrogen production and
264               !! that of silicon; depending upon the "intracellular" ratio
265               !! of Si:N, model diatoms will uptake nitrogen/silicon
266               !! differentially; this borrows from the diatom model of
267               !! Mongin et al. (2006)
268               !!----------------------------------------------------------
269               !!
270               if (jliebig .eq. 0) then
271                  !! multiplicative nutrient limitation
272                  fpdlim = fnld(ji,jj) * ffld(ji,jj)
273               elseif (jliebig .eq. 1) then
274                  !! Liebig Law (= most limiting) nutrient limitation
275                  fpdlim = min(fnld(ji,jj), ffld(ji,jj))
276               endif
277               !!
278          if (zphd(ji,jj).GT.rsmall .AND. zpds(ji,jj).GT.rsmall) then
279                  !! "intracellular" elemental ratios
280                  ! fsin(ji,jj)  = zpds(ji,jj) / (zphd(ji,jj) +              &
281                  !                               tiny(zphd(ji,jj)))
282                  ! fnsi         = zphd(ji,jj) / (zpds(ji,jj) +              &
283                  !                               tiny(zpds(ji,jj)))
284                  fsin(ji,jj) = 0.0
285                  IF( zphd(ji,jj) .GT. rsmall) fsin(ji,jj)  = zpds(ji,jj) /  &
286                                                              zphd(ji,jj)
287                  fnsi = 0.0
288                  IF( zpds(ji,jj) .GT. rsmall) fnsi  = zphd(ji,jj) /         &
289                                                       zpds(ji,jj)
290                  !! AXY (23/02/10): these next variables derive from
291                  !! Mongin et al. (2003)
292                  fsin1 = 3.0 * xsin0 !! = 0.6
293                  fnsi1 = 1.0 / fsin1 !! = 1.667
294                  fnsi2 = 1.0 / xsin0 !! = 5.0
295                  !!
296                  !! conditionalities based on ratios
297                  !! nitrogen (and iron and carbon)
298                  if (fsin(ji,jj).le.xsin0) then
299                     fprd(ji,jj)  = 0.0
300                     fsld2(ji,jj) = 0.0
301                  elseif (fsin(ji,jj).lt.fsin1) then
302                     fprd(ji,jj)  = xuif * ((fsin(ji,jj) - xsin0) /          &
303                                            (fsin(ji,jj) +                   &
304                                             tiny(fsin(ji,jj)))) *           &
305                                    (fjld(ji,jj) * fpdlim)
306                     fsld2(ji,jj) = xuif * ((fsin(ji,jj) - xsin0) /          &
307                                            (fsin(ji,jj) +                   &
308                                             tiny(fsin(ji,jj))))
309                  elseif (fsin(ji,jj).ge.fsin1) then
310                     fprd(ji,jj)  = (fjld(ji,jj) * fpdlim)
311                     fsld2(ji,jj) = 1.0
312                  endif
313                  !!
314                  !! silicon
315                  if (fsin(ji,jj).lt.fnsi1) then
316                     fprds(ji,jj) = (fjld(ji,jj) * fsld(ji,jj))
317                  elseif (fsin(ji,jj).lt.fnsi2) then
318                     fprds(ji,jj) = xuif * ((fnsi - xnsi0) /          &
319                                            (fnsi + tiny(fnsi))) *           &
320                                    (fjld(ji,jj) * fsld(ji,jj))
321                  else
322                     fprds(ji,jj) = 0.0
323                  endif     
324               else
325                  fsin(ji,jj)  = 0.0
326                  fnsi         = 0.0
327                  fprd(ji,jj)  = 0.0
328                  fsld2(ji,jj) = 0.0
329                  fprds(ji,jj) = 0.0
330               endif
331
332# if defined key_debug_medusa
333               !! report phytoplankton growth (including diatom silicon
334               !! submodel)
335               if (idf.eq.1.AND.idfval.eq.1) then
336                  IF (lwp) write (numout,*) '------------------------------'
337                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin(ji,jj)
338                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi
339                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2(ji,jj)
340                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn(ji,jj)
341                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd(ji,jj)
342                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds(ji,jj)
343               endif
344# endif
345            ENDIF
346         ENDDO
347      ENDDO
348
349      DO jj = 2,jpjm1
350         DO ji = 2,jpim1
351            if (tmask(ji,jj,jk) == 1) then
352               !!----------------------------------------------------------
353               !! Mixed layer primary production
354               !! this block calculates the amount of primary production
355               !! that occurs within the upper mixed layer; this allows the
356               !! separate diagnosis of "sub-surface" primary production; it
357               !! does assume that short-term variability in mixed layer
358               !! depth doesn't mess with things though
359               !!----------------------------------------------------------
360               !!
361               if (fdep1(ji,jj).le.hmld(ji,jj)) then
362                  !! this level is entirely in the mixed layer
363                  fq0 = 1.0
364               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
365                  !! this level is entirely below the mixed layer
366                  fq0 = 0.0
367               else
368                  !! this level straddles the mixed layer
369                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
370               endif
371               !!
372               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn(ji,jj) * zphn(ji,jj) * &
373                                                  fse3t(ji,jj,jk) * fq0)
374               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
375                                                  fse3t(ji,jj,jk) * fq0)
376          !! AXY (16/08/17)
377          fchl_ml(ji,jj) = fchl_ml(ji,jj) + ((zchn(ji,jj) + zchd(ji,jj)) * &
378                                             (fse3t(ji,jj,jk) * fq0) / hmld(ji,jj))
379            ENDIF
380         ENDDO
381      ENDDO
382
383      DO jj = 2,jpjm1
384         DO ji = 2,jpim1
385            if (tmask(ji,jj,jk) == 1) then
386               !!----------------------------------------------------------
387               !! Vertical Integral --
388               !!----------------------------------------------------------
389               !! vertical integral non-diatom phytoplankton
390               ftot_pn(ji,jj)  = ftot_pn(ji,jj)  + (zphn(ji,jj) *             &
391                                                    fse3t(ji,jj,jk))
392               !! vertical integral diatom phytoplankton
393               ftot_pd(ji,jj)  = ftot_pd(ji,jj)  + (zphd(ji,jj) *             &
394                                                    fse3t(ji,jj,jk))
395               !! vertical integral microzooplankton
396               ftot_zmi(ji,jj) = ftot_zmi(ji,jj) + (zzmi(ji,jj) *             &
397                                                    fse3t(ji,jj,jk))
398               !! vertical integral mesozooplankton
399               ftot_zme(ji,jj) = ftot_zme(ji,jj) + (zzme(ji,jj) *             &
400                                                    fse3t(ji,jj,jk))
401               !! vertical integral slow detritus, nitrogen
402               ftot_det(ji,jj) = ftot_det(ji,jj) + (zdet(ji,jj) *             &
403                                                    fse3t(ji,jj,jk))
404               !! vertical integral slow detritus, carbon
405               ftot_dtc(ji,jj) = ftot_dtc(ji,jj) + (zdtc(ji,jj) *             &
406                                                    fse3t(ji,jj,jk))
407            ENDIF
408         ENDDO
409      ENDDO
410
411      DO jj = 2,jpjm1
412         DO ji = 2,jpim1
413            if (tmask(ji,jj,jk) == 1) then
414               !!----------------------------------------------------------
415               !! More chlorophyll calculations
416               !!----------------------------------------------------------
417               !!
418               !! frn(ji,jj) = (xthetam / fthetan(ji,jj)) *                   &
419               !!              (fprn(ji,jj) / (fthetan(ji,jj) * xpar(ji,jj,jk)))
420               !! frd(ji,jj) = (xthetam / fthetad(ji,jj)) *                   &
421               !!              (fprd(ji,jj) / (fthetad(ji,jj) * xpar(ji,jj,jk)))
422               frn(ji,jj) = (xthetam * fchn(ji,jj) * fnln(ji,jj) *            &
423                             ffln2(ji,jj)) / (fthetan(ji,jj) +                &
424                                             tiny(fthetan(ji,jj)))
425               !! AXY (12/05/09): there's potentially a problem here; fsld,
426               !!   silicic acid limitation, is used in the following line
427               !!   to regulate chlorophyll growth in a manner that is
428               !!   inconsistent with its use in the regulation of biomass
429               !!   growth; the Mongin term term used in growth is more
430               !!   complex than the simple multiplicative function used
431               !!   below
432               !! frd(ji,jj) = (xthetam * fchd(ji,jj) * fnld(ji,jj) *        &
433               !!               ffld(ji,jj) * fsld(ji,jj)) /                 &
434               !!               (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
435               !! AXY (12/05/09): this replacement line uses the new
436               !!   variable, fsld2, to regulate chlorophyll growth
437               frd(ji,jj) = (xthetamd * fchd(ji,jj) * fnld(ji,jj) *          &
438                             ffld(ji,jj) * fsld2(ji,jj)) /                   &
439                             (fthetad(ji,jj) + tiny(fthetad(ji,jj)))
440
441# if defined key_debug_medusa
442               !! report chlorophyll calculations
443               if (idf.eq.1.AND.idfval.eq.1) then
444                  IF (lwp) write (numout,*) '------------------------------'
445                  IF (lwp) write (numout,*) 'fthetan(',jk,') = ', fthetan(ji,jj)
446                  IF (lwp) write (numout,*) 'fthetad(',jk,') = ', fthetad(ji,jj)
447                  IF (lwp) write (numout,*) 'frn(',jk,')     = ', frn(ji,jj)
448                  IF (lwp) write (numout,*) 'frd(',jk,')     = ', frd(ji,jj)
449               endif
450# endif
451
452            ENDIF
453         ENDDO
454      ENDDO
455
456   END SUBROUTINE phytoplankton
457
458#else
459   !!======================================================================
460   !!  Dummy module :                                   No MEDUSA bio-model
461   !!======================================================================
462CONTAINS
463   SUBROUTINE phytoplankton( )                    ! Empty routine
464      WRITE(*,*) 'phytoplankton: You should not have seen this print! error?'
465   END SUBROUTINE phytoplankton
466#endif 
467
468   !!======================================================================
469END MODULE phytoplankton_mod
Note: See TracBrowser for help on using the repository browser.