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/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/UKMO/dev_r5518_medusa_chg_trc_bio_medusa/NEMOGCM/NEMO/TOP_SRC/MEDUSA/phytoplankton.F90 @ 7975

Last change on this file since 7975 was 7975, checked in by marc, 7 years ago

Removed plankton processes from trcbio_medusa.F90 into extra routines

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   !!----------------------------------------------------------------------
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: fald, faln, fchd, fchd1, fchn, fchn1, &
34                                   fdep1, ffld, ffln2, fjld, fjln,       &
35                                   fjlim_pd, fjlim_pn, fjln,             &
36                                   fnld, fnln, fnsi,                     &
37                                   fpdlim, fpnlim, 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                                   xvpdT, xvpnT,                         &
44                                   zchd, zchn, zdet, zdin, zdtc,         &
45                                   zfer, zpds, zphd, zphn, zsil,         &
46                                   zzme, zzmi
47      USE dom_oce,           ONLY: e3t_0, e3t_n, gdepw_0, gdepw_n, tmask
48      USE in_out_manager,    ONLY: lwp, numout
49      USE oce,               ONLY: tsn
50      USE par_kind,          ONLY: wp
51      USE par_oce,           ONLY: jp_tem, jpim1, jpjm1
52      USE phycst,            ONLY: rsmall
53      USE sms_medusa,        ONLY: jliebig, jphy, jq10,                  &
54                                   xald, xaln, xfld, xfln,               &
55                                   xnld, xnln, xnsi0, xpar,              &
56                                   xsin0, xsld, xthetam, xthetamd, xuif, &
57                                   xvpd, xvpn, xxi
58      USE zdfmxl,            ONLY: hmld
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)                     ::    fsin1,fnsi1,fnsi2
69      REAL(wp) ::    fq0
70
71      DO jj = 2,jpjm1
72         DO ji = 2,jpim1
73            !! OPEN wet point IF..THEN loop
74            if (tmask(ji,jj,jk) == 1) then
75               !!----------------------------------------------------------
76               !! Chlorophyll calculations
77               !!----------------------------------------------------------
78               !!
79               !! non-diatoms
80          if (zphn(ji,jj).GT.rsmall) then
81                  fthetan(ji,jj) = max(tiny(zchn(ji,jj)),                    &
82                                       (zchn(ji,jj) * xxi) /                 &
83                                       (zphn(ji,jj) + tiny(zphn(ji,jj))))
84                  faln(ji,jj)    = xaln * fthetan(ji,jj)
85               else
86                  fthetan(ji,jj) = 0.
87                  faln(ji,jj)    = 0.
88               endif
89               !!
90               !! diatoms
91          if (zphd(ji,jj).GT.rsmall) then
92                  fthetad(ji,jj) = max(tiny(zchd(ji,jj)),                   &
93                                       (zchd(ji,jj) * xxi) /                &
94                                       (zphd(ji,jj) + tiny(zphd(ji,jj))))
95                  fald(ji,jj)    = xald * fthetad(ji,jj)
96               else
97                  fthetad(ji,jj) = 0.
98                  fald(ji,jj)    = 0.
99               endif
100
101# if defined key_debug_medusa
102               !! report biological calculations
103               if (idf.eq.1.AND.idfval.eq.1) then
104                  IF (lwp) write (numout,*) '------------------------------'
105                  IF (lwp) write (numout,*) 'faln(',jk,') = ', faln(ji,jj)
106                  IF (lwp) write (numout,*) 'fald(',jk,') = ', fald(ji,jj)
107               endif
108# endif
109            ENDIF
110         ENDDO
111      ENDDO
112
113      DO jj = 2,jpjm1
114         DO ji = 2,jpim1
115            if (tmask(ji,jj,1) == 1) then
116               !!----------------------------------------------------------
117               !! Phytoplankton light limitation
118               !!----------------------------------------------------------
119               !!
120               !! It is assumed xpar is the depth-averaged (vertical layer) PAR
121               !! Light limitation (check self-shading) in W/m2
122               !!
123               !! Note that there is no temperature dependence in phytoplankton
124               !! growth rate or any other function.
125               !! In calculation of Chl/Phy ratio tiny(phyto) is introduced to
126               !! avoid NaNs in case of Phy==0. 
127               !!
128               !! fthetad and fthetan are Chl:C ratio (gChl/gC) in diat and
129               !! non-diat:
130               !! for 1:1 Chl:P ratio (mgChl/mmolN) theta=0.012
131               !!
132               !! AXY (16/07/09)
133               !! temperature for new Eppley style phytoplankton growth
134               fun_T(ji,jj)   = 1.066**(1.0 * tsn(ji,jj,jk,jp_tem))
135               !! AXY (16/05/11): add in new Q10 (1.5, not 2.0) for
136               !!                 phytoplankton growth; remin. unaffected
137               fun_Q10(ji,jj) = jq10**((tsn(ji,jj,jk,jp_tem) - 0.0) / 10.0)
138               if (jphy.eq.1) then
139                  xvpnT(ji,jj) = xvpn * fun_T(ji,jj)
140                  xvpdT(ji,jj) = xvpd * fun_T(ji,jj)
141               elseif (jphy.eq.2) then
142                  xvpnT(ji,jj) = xvpn * fun_Q10(ji,jj)
143                  xvpdT(ji,jj) = xvpd * fun_Q10(ji,jj)
144               else
145                  xvpnT(ji,jj) = xvpn
146                  xvpdT(ji,jj) = xvpd
147               endif
148            ENDIF
149         ENDDO
150      ENDDO
151
152      DO jj = 2,jpjm1
153         DO ji = 2,jpim1
154            if (tmask(ji,jj,1) == 1) then           
155               !!
156               !! non-diatoms
157               fchn1(ji,jj)   = (xvpnT(ji,jj) * xvpnT(ji,jj)) +              &
158                                 (faln(ji,jj) * faln(ji,jj) *                &
159                                  xpar(ji,jj,jk) * xpar(ji,jj,jk))
160               if (fchn1(ji,jj).GT.rsmall) then
161                  fchn(ji,jj)    = xvpnT(ji,jj) / (sqrt(fchn1(ji,jj)) +      &
162                                                   tiny(fchn1(ji,jj)))
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(ji,jj)
169            ENDIF
170         ENDDO
171      ENDDO
172
173      DO jj = 2,jpjm1
174         DO ji = 2,jpim1
175            if (tmask(ji,jj,1) == 1) then
176               !!
177               !! diatoms
178               fchd1(ji,jj)   = (xvpdT(ji,jj) * xvpdT(ji,jj)) +              &
179                                (fald(ji,jj) * fald(ji,jj) *                 &
180                                 xpar(ji,jj,jk) * xpar(ji,jj,jk))
181               if (fchd1(ji,jj).GT.rsmall) then
182                  fchd(ji,jj)    = xvpdT(ji,jj) / (sqrt(fchd1(ji,jj)) +      &
183                                                   tiny(fchd1(ji,jj)))
184               else
185                  fchd(ji,jj)    = 0.
186               endif
187               !! diatom J term
188               fjld(ji,jj)    = fchd(ji,jj) * fald(ji,jj) * xpar(ji,jj,jk)
189               fjlim_pd(ji,jj) = fjld(ji,jj) / xvpdT(ji,jj)
190     
191# if defined key_debug_medusa
192               !! report phytoplankton light limitation
193               if (idf.eq.1.AND.idfval.eq.1) then
194                  IF (lwp) write (numout,*) '------------------------------'
195                  IF (lwp) write (numout,*) 'fchn(',jk,') = ', fchn(ji,jj)
196                  IF (lwp) write (numout,*) 'fchd(',jk,') = ', fchd(ji,jj)
197                  IF (lwp) write (numout,*) 'fjln(',jk,') = ', fjln(ji,jj)
198                  IF (lwp) write (numout,*) 'fjld(',jk,') = ', fjld(ji,jj)
199               endif
200# endif
201            ENDIF
202         ENDDO
203      ENDDO
204
205      DO jj = 2,jpjm1
206         DO ji = 2,jpim1
207            if (tmask(ji,jj,1) == 1) then
208               !!----------------------------------------------------------
209               !! Phytoplankton nutrient limitation
210               !!----------------------------------------------------------
211               !!
212               !! non-diatoms (N, Fe).
213               !! non-diatom Qn term
214               fnln(ji,jj)  = zdin(ji,jj) / (zdin(ji,jj) + xnln)
215               !! non-diatom Qf term
216               ffln2(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfln)
217               !!
218               !! diatoms (N, Si, Fe).
219               !! diatom Qn term
220               fnld(ji,jj) = zdin(ji,jj) / (zdin(ji,jj) + xnld)
221               !! diatom Qs term
222               fsld(ji,jj) = zsil(ji,jj) / (zsil(ji,jj) + xsld)
223               !! diatom Qf term
224               ffld(ji,jj) = zfer(ji,jj) / (zfer(ji,jj) + xfld)
225
226# if defined key_debug_medusa
227               !! report phytoplankton nutrient limitation
228               if (idf.eq.1.AND.idfval.eq.1) then
229                  IF (lwp) write (numout,*) '------------------------------'
230                  IF (lwp) write (numout,*) 'fnln(',jk,') = ', fnln(ji,jj)
231                  IF (lwp) write (numout,*) 'fnld(',jk,') = ', fnld(ji,jj)
232                  IF (lwp) write (numout,*) 'ffln2(',jk,') = ', ffln2(ji,jj)
233                  IF (lwp) write (numout,*) 'ffld(',jk,') = ', ffld(ji,jj)
234                  IF (lwp) write (numout,*) 'fsld(',jk,') = ', fsld(ji,jj)
235               endif
236# endif
237            ENDIF
238         ENDDO
239      ENDDO
240
241      DO jj = 2,jpjm1
242         DO ji = 2,jpim1
243            if (tmask(ji,jj,1) == 1) then
244               !!----------------------------------------------------------
245               !! Primary production (non-diatoms)
246               !! (note: still needs multiplying by phytoplankton
247               !! concentration)
248               !!----------------------------------------------------------
249               !!
250               if (jliebig .eq. 0) then
251                  !! multiplicative nutrient limitation
252                  fpnlim(ji,jj) = fnln(ji,jj) * ffln2(ji,jj)
253               elseif (jliebig .eq. 1) then
254                  !! Liebig Law (= most limiting) nutrient limitation
255                  fpnlim(ji,jj) = min(fnln(ji,jj), ffln2(ji,jj))
256               endif
257               fprn(ji,jj) = fjln(ji,jj) * fpnlim(ji,jj)
258            ENDIF
259         ENDDO
260      ENDDO
261
262      DO jj = 2,jpjm1
263         DO ji = 2,jpim1
264            if (tmask(ji,jj,1) == 1) then
265               !!----------------------------------------------------------
266               !! Primary production (diatoms)
267               !! (note: still needs multiplying by phytoplankton
268               !! concentration)
269               !!
270               !! Production here is split between nitrogen production and
271               !! that of silicon; depending upon the "intracellular" ratio
272               !! of Si:N, model diatoms will uptake nitrogen/silicon
273               !! differentially; this borrows from the diatom model of
274               !! Mongin et al. (2006)
275               !!----------------------------------------------------------
276               !!
277               if (jliebig .eq. 0) then
278                  !! multiplicative nutrient limitation
279                  fpdlim(ji,jj) = fnld(ji,jj) * ffld(ji,jj)
280               elseif (jliebig .eq. 1) then
281                  !! Liebig Law (= most limiting) nutrient limitation
282                  fpdlim(ji,jj) = min(fnld(ji,jj), ffld(ji,jj))
283               endif
284               !!
285          if (zphd(ji,jj).GT.rsmall .AND. zpds(ji,jj).GT.rsmall) then
286                  !! "intracellular" elemental ratios
287                  ! fsin(ji,jj)  = zpds(ji,jj) / (zphd(ji,jj) +              &
288                  !                               tiny(zphd(ji,jj)))
289                  ! fnsi(ji,jj)  = zphd(ji,jj) / (zpds(ji,jj) +              &
290                  !                               tiny(zpds(ji,jj)))
291                  fsin(ji,jj) = 0.0
292                  IF( zphd(ji,jj) .GT. rsmall) fsin(ji,jj)  = zpds(ji,jj) /  &
293                                                              zphd(ji,jj)
294                  fnsi(ji,jj) = 0.0
295                  IF( zpds(ji,jj) .GT. rsmall) fnsi(ji,jj)  = zphd(ji,jj) /  &
296                                                              zpds(ji,jj)
297                  !! AXY (23/02/10): these next variables derive from
298                  !! Mongin et al. (2003)
299                  fsin1 = 3.0 * xsin0 !! = 0.6
300                  fnsi1 = 1.0 / fsin1 !! = 1.667
301                  fnsi2 = 1.0 / xsin0 !! = 5.0
302                  !!
303                  !! conditionalities based on ratios
304                  !! nitrogen (and iron and carbon)
305                  if (fsin(ji,jj).le.xsin0) then
306                     fprd(ji,jj)  = 0.0
307                     fsld2(ji,jj) = 0.0
308                  elseif (fsin(ji,jj).lt.fsin1) then
309                     fprd(ji,jj)  = xuif * ((fsin(ji,jj) - xsin0) /          &
310                                            (fsin(ji,jj) +                   &
311                                             tiny(fsin(ji,jj)))) *           &
312                                    (fjld(ji,jj) * fpdlim(ji,jj))
313                     fsld2(ji,jj) = xuif * ((fsin(ji,jj) - xsin0) /          &
314                                            (fsin(ji,jj) +                   &
315                                             tiny(fsin(ji,jj))))
316                  elseif (fsin(ji,jj).ge.fsin1) then
317                     fprd(ji,jj)  = (fjld(ji,jj) * fpdlim(ji,jj))
318                     fsld2(ji,jj) = 1.0
319                  endif
320                  !!
321                  !! silicon
322                  if (fsin(ji,jj).lt.fnsi1) then
323                     fprds(ji,jj) = (fjld(ji,jj) * fsld(ji,jj))
324                  elseif (fsin(ji,jj).lt.fnsi2) then
325                     fprds(ji,jj) = xuif * ((fnsi(ji,jj) - xnsi0) /          &
326                                            (fnsi(ji,jj) +                   &
327                                             tiny(fnsi(ji,jj)))) *           &
328                                    (fjld(ji,jj) * fsld(ji,jj))
329                  else
330                     fprds(ji,jj) = 0.0
331                  endif     
332               else
333                  fsin(ji,jj)  = 0.0
334                  fnsi(ji,jj)  = 0.0
335                  fprd(ji,jj)  = 0.0
336                  fsld2(ji,jj) = 0.0
337                  fprds(ji,jj) = 0.0
338               endif
339
340# if defined key_debug_medusa
341               !! report phytoplankton growth (including diatom silicon
342               !! submodel)
343               if (idf.eq.1.AND.idfval.eq.1) then
344                  IF (lwp) write (numout,*) '------------------------------'
345                  IF (lwp) write (numout,*) 'fsin(',jk,')   = ', fsin(ji,jj)
346                  IF (lwp) write (numout,*) 'fnsi(',jk,')   = ', fnsi(ji,jj)
347                  IF (lwp) write (numout,*) 'fsld2(',jk,')  = ', fsld2(ji,jj)
348                  IF (lwp) write (numout,*) 'fprn(',jk,')   = ', fprn(ji,jj)
349                  IF (lwp) write (numout,*) 'fprd(',jk,')   = ', fprd(ji,jj)
350                  IF (lwp) write (numout,*) 'fprds(',jk,')  = ', fprds(ji,jj)
351               endif
352# endif
353            ENDIF
354         ENDDO
355      ENDDO
356
357      DO jj = 2,jpjm1
358         DO ji = 2,jpim1
359            if (tmask(ji,jj,1) == 1) then
360               !!----------------------------------------------------------
361               !! Mixed layer primary production
362               !! this block calculates the amount of primary production
363               !! that occurs within the upper mixed layer; this allows the
364               !! separate diagnosis of "sub-surface" primary production; it
365               !! does assume that short-term variability in mixed layer
366               !! depth doesn't mess with things though
367               !!----------------------------------------------------------
368               !!
369               if (fdep1(ji,jj).le.hmld(ji,jj)) then
370                  !! this level is entirely in the mixed layer
371                  fq0 = 1.0
372               elseif (fsdepw(ji,jj,jk).ge.hmld(ji,jj)) then
373                  !! this level is entirely below the mixed layer
374                  fq0 = 0.0
375               else
376                  !! this level straddles the mixed layer
377                  fq0 = (hmld(ji,jj) - fsdepw(ji,jj,jk)) / fse3t(ji,jj,jk)
378               endif
379               !!
380               fprn_ml(ji,jj) = fprn_ml(ji,jj) + (fprn(ji,jj) * zphn(ji,jj) * &
381                                                  fse3t(ji,jj,jk) * fq0)
382               fprd_ml(ji,jj) = fprd_ml(ji,jj) + (fprd(ji,jj) * zphd(ji,jj) * &
383                                                  fse3t(ji,jj,jk) * fq0)
384            ENDIF
385         ENDDO
386      ENDDO
387
388      DO jj = 2,jpjm1
389         DO ji = 2,jpim1
390            if (tmask(ji,jj,1) == 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,1) == 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.