source: branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90 @ 9309

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

Apply changes developed under Met Office GMED ticket number 351 in development
branch branches/NERC/dev_r5518_GO6_ScalingCoupledChl.

The command issued to perform the merge is:

svn merge -r 8590:9053 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/NERC/dev_r5518_GO6_ScalingCoupledChl

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