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

Last change on this file since 10149 was 10020, checked in by marc, 2 years ago

GMED ticket 406. CPP key fixes.

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