source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90 @ 8441

Last change on this file since 8441 was 8441, checked in by frrh, 3 years ago

Commit changes relating to Met Office GMED ticket 339 for the modularisation of
of trcbio_medusa.F90.

Branch http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
from revisions 8394 to 8423 inclusive refer.

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