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.
p5zprod.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/TOP/PISCES/P4Z/p5zprod.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

Add Mixed Precision support by Oriol Tintó
  • Property svn:keywords set to Id
File size: 31.8 KB
Line 
1MODULE p5zprod
2   !!======================================================================
3   !!                         ***  MODULE p5zprod  ***
4   !! TOP :  Growth Rate of the two phytoplanktons groups
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-05  (O. Aumont, C. Ethe) New parameterization of light limitation
9   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
10   !!----------------------------------------------------------------------
11   !!   p5z_prod       :   Compute the growth Rate of the two phytoplanktons groups
12   !!   p5z_prod_init  :   Initialization of the parameters for growth
13   !!   p5z_prod_alloc :   Allocate variables for growth
14   !!----------------------------------------------------------------------
15   USE oce_trc         !  shared variables between ocean and passive tracers
16   USE trc             !  passive tracers common variables
17   USE sms_pisces      !  PISCES Source Minus Sink variables
18   USE p4zlim
19   USE p5zlim          !  Co-limitations of differents nutrients
20   USE prtctl          !  print control for debugging
21   USE iom             !  I/O manager
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   p5z_prod         ! called in p5zbio.F90
27   PUBLIC   p5z_prod_init    ! called in trcsms_pisces.F90
28   PUBLIC   p5z_prod_alloc
29
30   !! * Shared module variables
31   REAL(wp), PUBLIC ::  pislopen        !:
32   REAL(wp), PUBLIC ::  pislopep        !:
33   REAL(wp), PUBLIC ::  pisloped        !:
34   REAL(wp), PUBLIC ::  xadap           !:
35   REAL(wp), PUBLIC ::  excretn         !:
36   REAL(wp), PUBLIC ::  excretp         !:
37   REAL(wp), PUBLIC ::  excretd         !:
38   REAL(wp), PUBLIC ::  bresp           !:
39   REAL(wp), PUBLIC ::  thetanpm        !:
40   REAL(wp), PUBLIC ::  thetannm        !:
41   REAL(wp), PUBLIC ::  thetandm        !:
42   REAL(wp), PUBLIC ::  chlcmin         !:
43   REAL(wp), PUBLIC ::  grosip          !:
44
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   zdaylen
46   
47   REAL(wp) :: r1_rday                !: 1 / rday
48   REAL(wp) :: texcretn               !: 1 - excret
49   REAL(wp) :: texcretp               !: 1 - excretp
50   REAL(wp) :: texcretd               !: 1 - excret2       
51
52   !! * Substitutions
53#  include "do_loop_substitute.h90"
54#  include "domzgr_substitute.h90"
55#  include "single_precision_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
58   !! $Id$
59   !! Software governed by the CeCILL license (see ./LICENSE)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE p5z_prod( kt , knt, Kbb, Kmm, Krhs )
64      !!---------------------------------------------------------------------
65      !!                     ***  ROUTINE p5z_prod  ***
66      !!
67      !! ** Purpose :   Compute the phytoplankton production depending on
68      !!              light, temperature and nutrient availability
69      !!
70      !! ** Method  : - ???
71      !!---------------------------------------------------------------------
72      !
73      INTEGER, INTENT(in) :: kt, knt
74      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs      ! time level indices
75      !
76      INTEGER  ::   ji, jj, jk
77      REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
78      REAL(wp) ::   zration, zratiop, zratiof, zmax, zmax2, zsilim, ztn, zadap
79      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
80      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
81      REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
82      REAL(wp) ::   zpislopen, zpislopep, zpisloped, thetannm_n, thetandm_n, thetanpm_n
83      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
84      REAL(wp) ::   zfact, zrfact2
85      CHARACTER (len=25) :: charout
86      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat, zstrn
87      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
88      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
89      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprbio, zprpic, zprdia, zysopt
90      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprchln, zprchlp, zprchld
91      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprorcan, zprorcap, zprorcad 
92      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprofed, zprofep, zprofen
93      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpronewn, zpronewp, zpronewd
94      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zproregn, zproregp, zproregd
95      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpropo4n, zpropo4p, zpropo4d
96      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprodopn, zprodopp, zprodopd
97      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrespn, zrespp, zrespd
98      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcroissn, zcroissp, zcroissd
99      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxl_fac, zmxl_chl
100      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpligprod1, zpligprod2
101      !!---------------------------------------------------------------------
102      !
103      IF( ln_timing )   CALL timing_start('p5z_prod')
104      !
105      zprorcan(:,:,:) = 0._wp ; zprorcap(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp
106      zcroissn(:,:,:) = 0._wp ; zcroissp(:,:,:) = 0._wp ; zcroissd(:,:,:) = 0._wp
107      zprofed (:,:,:) = 0._wp ; zprofep (:,:,:) = 0._wp ; zprofen (:,:,:) = 0._wp
108      zpronewn(:,:,:) = 0._wp ; zpronewp(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp
109      zproregn(:,:,:) = 0._wp ; zproregp(:,:,:) = 0._wp ; zproregd(:,:,:) = 0._wp 
110      zpropo4n(:,:,:) = 0._wp ; zpropo4p(:,:,:) = 0._wp ; zpropo4d(:,:,:) = 0._wp
111      zprdia  (:,:,:) = 0._wp ; zprpic  (:,:,:) = 0._wp ; zprbio  (:,:,:) = 0._wp
112      zprodopn(:,:,:) = 0._wp ; zprodopp(:,:,:) = 0._wp ; zprodopd(:,:,:) = 0._wp
113      zysopt  (:,:,:) = 0._wp 
114      zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
115
116      ! Computation of the optimal production
117      zprnut (:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
118      zprmaxn(:,:,:) = ( 0.65_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc(:,:,:)
119      zprmaxp(:,:,:) = 0.5 / 0.65 * zprmaxn(:,:,:) 
120      zprmaxd(:,:,:) = zprmaxn(:,:,:) 
121
122      ! compute the day length depending on latitude and the day
123      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
124      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
125
126      ! day length in hours
127      zstrn(:,:) = 0.
128      DO_2D( 1, 1, 1, 1 )
129         zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
130         zargu = MAX( -1., MIN(  1., zargu ) )
131         zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
132      END_2D
133
134         ! Impact of the day duration on phytoplankton growth
135      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
136         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
137            zval = MAX( 1., zstrn(ji,jj) )
138            IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
139               zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
140            ENDIF
141            zmxl_chl(ji,jj,jk) = zval / 24.
142            zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval )
143         ENDIF
144      END_3D
145
146      zprbio(:,:,:) = zprmaxn(:,:,:) * zmxl_fac(:,:,:)
147      zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
148      zprpic(:,:,:) = zprmaxp(:,:,:) * zmxl_fac(:,:,:)
149
150
151      ! Maximum light intensity
152      zdaylen(:,:) = MAX(1., zstrn(:,:)) / 24.
153      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
154
155      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
156         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
157            ! Computation of the P-I slope for nanos and diatoms
158            ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
159            zadap       = xadap * ztn / ( 2.+ ztn )
160            !
161            zpislopeadn(ji,jj,jk) = pislopen * tr(ji,jj,jk,jpnch,Kbb)    &
162            &                       /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn)
163            zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.25 * epico(ji,jj,jk) ) )   &
164            &                       * tr(ji,jj,jk,jppch,Kbb) /( tr(ji,jj,jk,jppic,Kbb) * 12. + rtrn)
165            zpislopeadd(ji,jj,jk) = pisloped * tr(ji,jj,jk,jpdch,Kbb)    &
166               &                    /( tr(ji,jj,jk,jpdia,Kbb) * 12. + rtrn)
167            !
168            zpislopen = zpislopeadn(ji,jj,jk) / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
169            zpislopep = zpislopeadp(ji,jj,jk) / ( zprpic(ji,jj,jk) * rday * xlimpic(ji,jj,jk) + rtrn )
170            zpisloped = zpislopeadd(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
171
172            ! Computation of production function for Carbon
173            !  ---------------------------------------------
174            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
175            zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
176            zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
177
178            ! Computation of production function for Chlorophyll
179            !  -------------------------------------------------
180            zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
181            zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
182            zpislopep = zpislopep * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
183            zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) )  )
184            zprchlp(ji,jj,jk) = zprmaxp(ji,jj,jk) * ( 1.- EXP( -zpislopep * epicom(ji,jj,jk) )  )
185            zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) )  )
186         ENDIF
187      END_3D
188
189      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
190
191          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
192            !    Si/C of diatoms
193            !    ------------------------
194            !    Si/C increases with iron stress and silicate availability
195            !    Si/C is arbitrariliy increased for very high Si concentrations
196            !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
197            zlim  = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi1 )
198            zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
199            zsilfac = 3.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
200            zsiborn = tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb)
201            IF (gphit(ji,jj) < -30 ) THEN
202              zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
203            ELSE
204              zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
205            ENDIF
206            zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
207        ENDIF
208      END_3D
209
210      !  Sea-ice effect on production                                                                               
211      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
212         zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
213         zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
214         zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
215         zprnut(ji,jj,jk)  = zprnut(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
216      END_3D
217
218      ! Computation of the various production terms of nanophytoplankton
219      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
220         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
221            !  production terms for nanophyto.
222            zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2
223            !
224            zration = tr(ji,jj,jk,jpnph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn )
225            zratiop = tr(ji,jj,jk,jppph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn )
226            zratiof = tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn )
227            zprnutmax = zprnut(ji,jj,jk) * fvnuptk(ji,jj,jk) / rno3 * tr(ji,jj,jk,jpphy,Kbb) * rfact2
228            ! Uptake of nitrogen
229            zrat = MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) ) 
230            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
231            zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpnmin(ji,jj,jk) )   &
232            &          / ( xqpnmax(ji,jj,jk) - xqpnmin(ji,jj,jk) + rtrn ), xlimnfe(ji,jj,jk) ) )
233            zpronewn(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xnanono3(ji,jj,jk)
234            zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
235            ! Uptake of phosphorus
236            zrat = MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
237            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
238            zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk)
239            zpropo4n(ji,jj,jk) = zpropmax * xnanopo4(ji,jj,jk)
240            zprodopn(ji,jj,jk) = zpropmax * xnanodop(ji,jj,jk)
241            ! Uptake of iron
242            zrat = MIN( 1., zratiof / qfnmax )
243            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
244            zprofmax = zprnutmax * qfnmax * zmax
245            zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk) * ( 3. - 2.4 * xlimnfe(ji,jj,jk)    &
246            &          / ( xlimnfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
247            &          + xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )
248         ENDIF
249      END_3D
250
251      ! Computation of the various production terms of picophytoplankton
252      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
253         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
254            !  production terms for picophyto.
255            zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * tr(ji,jj,jk,jppic,Kbb) * rfact2
256            !
257            zration = tr(ji,jj,jk,jpnpi,Kbb) / ( tr(ji,jj,jk,jppic,Kbb) + rtrn )
258            zratiop = tr(ji,jj,jk,jpppi,Kbb) / ( tr(ji,jj,jk,jppic,Kbb) + rtrn )
259            zratiof = tr(ji,jj,jk,jppfe,Kbb) / ( tr(ji,jj,jk,jppic,Kbb) + rtrn )
260            zprnutmax = zprnut(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * tr(ji,jj,jk,jppic,Kbb) * rfact2
261            ! Uptake of nitrogen
262            zrat = MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
263            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
264            zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqppmin(ji,jj,jk) )   &
265            &          / ( xqppmax(ji,jj,jk) - xqppmin(ji,jj,jk) + rtrn ), xlimpfe(ji,jj,jk) ) )
266            zpronewp(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xpicono3(ji,jj,jk) 
267            zproregp(ji,jj,jk) = zpronmax * xpiconh4(ji,jj,jk)
268            ! Uptake of phosphorus
269            zrat = MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
270            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
271            zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk)
272            zpropo4p(ji,jj,jk) = zpropmax * xpicopo4(ji,jj,jk)
273            zprodopp(ji,jj,jk) = zpropmax * xpicodop(ji,jj,jk)
274            ! Uptake of iron
275            zrat = MIN( 1., zratiof / qfpmax )
276            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
277            zprofmax = zprnutmax * qfpmax * zmax
278            zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk) * ( 3. - 2.4 * xlimpfe(ji,jj,jk)   &
279            &          / ( xlimpfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
280            &          + xpicono3(ji,jj,jk) + xpiconh4(ji,jj,jk) ) * (1. - xpicofer(ji,jj,jk) ) )
281         ENDIF
282      END_3D
283
284      ! Computation of the various production terms of diatoms
285      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
286         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
287            !  production terms for diatomees
288            zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * tr(ji,jj,jk,jpdia,Kbb) * rfact2
289            ! Computation of the respiration term according to pahlow
290            ! & oschlies (2013)
291            !
292            zration = tr(ji,jj,jk,jpndi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
293            zratiop = tr(ji,jj,jk,jppdi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
294            zratiof = tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
295            zprnutmax = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * tr(ji,jj,jk,jpdia,Kbb) * rfact2
296            ! Uptake of nitrogen
297            zrat = MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
298            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
299            zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpdmin(ji,jj,jk) )   &
300            &          / ( xqpdmax(ji,jj,jk) - xqpdmin(ji,jj,jk) + rtrn ), xlimdfe(ji,jj,jk) ) )
301            zpronewd(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xdiatno3(ji,jj,jk)
302            zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
303            ! Uptake of phosphorus
304            zrat = MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
305            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
306            zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
307            zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
308            zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
309            ! Uptake of iron
310            zrat = MIN( 1., zratiof / qfdmax )
311            zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05))
312            zprofmax = zprnutmax * qfdmax * zmax
313            zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk) * ( 3. - 2.4 * xlimdfe(ji,jj,jk)     &
314            &          / ( xlimdfe(ji,jj,jk) + 0.2 ) ) * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
315            &          + xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )
316         ENDIF
317      END_3D
318
319      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
320         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
321               !  production terms for nanophyto. ( chlorophyll )
322            znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
323            zprod = rday * (zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
324            thetannm_n   = MIN ( thetannm, ( thetannm / (1. - 1.14 / 43.4 *ts(ji,jj,jk,jp_tem,Kmm)))   &
325            &               * (1. - 1.14 / 43.4 * 20.))
326            zprochln = thetannm_n * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
327            zprochln = MAX(zprochln, chlcmin * 12. * zprorcan (ji,jj,jk) )
328               !  production terms for picophyto. ( chlorophyll )
329            zpicotot = epicom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
330            zprod = rday * (zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)) * zprchlp(ji,jj,jk) * xlimpic(ji,jj,jk)
331            thetanpm_n   = MIN ( thetanpm, ( thetanpm / (1. - 1.14 / 43.4 *ts(ji,jj,jk,jp_tem,Kmm)))   &
332            &               * (1. - 1.14 / 43.4 * 20.))
333            zprochlp = thetanpm_n * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
334            zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
335            !  production terms for diatomees ( chlorophyll )
336            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
337            zprod = rday * (zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
338            thetandm_n   = MIN ( thetandm, ( thetandm / (1. - 1.14 / 43.4 *ts(ji,jj,jk,jp_tem,Kmm)))   &
339            &               * (1. - 1.14 / 43.4 * 20.))
340            zprochld = thetandm_n * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
341            zprochld = MAX(zprochld, chlcmin * 12. * zprorcad(ji,jj,jk) )
342            !   Update the arrays TRA which contain the Chla sources and sinks
343            tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) + zprochln * texcretn
344            tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) + zprochld * texcretd
345            tr(ji,jj,jk,jppch,Krhs) = tr(ji,jj,jk,jppch,Krhs) + zprochlp * texcretp
346         ENDIF
347      END_3D
348
349      !   Update the arrays TRA which contain the biological sources and sinks
350      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
351        zprontot = zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)
352        zproptot = zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)
353        zprodtot = zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)
354        zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
355        &          + excretp * zprorcap(ji,jj,jk)
356        tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zpropo4n(ji,jj,jk) - zpropo4d(ji,jj,jk)  &
357        &                     - zpropo4p(ji,jj,jk)
358        tr(ji,jj,jk,jpno3,Krhs) = tr(ji,jj,jk,jpno3,Krhs) - zpronewn(ji,jj,jk) - zpronewd(ji,jj,jk)  &
359        &                     - zpronewp(ji,jj,jk)
360        tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) - zproregn(ji,jj,jk) - zproregd(ji,jj,jk)  &
361        &                     - zproregp(ji,jj,jk)
362        tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) + zprorcan(ji,jj,jk) * texcretn    &
363           &                  - zpsino3 * zpronewn(ji,jj,jk) - zpsinh4 * zproregn(ji,jj,jk)   &
364           &                  - zrespn(ji,jj,jk) 
365        zcroissn(ji,jj,jk) = tr(ji,jj,jk,jpphy,Krhs) / rfact2/ (tr(ji,jj,jk,jpphy,Kbb) + rtrn)
366        tr(ji,jj,jk,jpnph,Krhs) = tr(ji,jj,jk,jpnph,Krhs) + zprontot * texcretn
367        tr(ji,jj,jk,jppph,Krhs) = tr(ji,jj,jk,jppph,Krhs) + zpropo4n(ji,jj,jk) * texcretn   &
368        &                     + zprodopn(ji,jj,jk) * texcretn
369        tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) + zprofen(ji,jj,jk) * texcretn
370        tr(ji,jj,jk,jppic,Krhs) = tr(ji,jj,jk,jppic,Krhs) + zprorcap(ji,jj,jk) * texcretp     &
371           &                  - zpsino3 * zpronewp(ji,jj,jk) - zpsinh4 * zproregp(ji,jj,jk)   &
372           &                  - zrespp(ji,jj,jk) 
373        zcroissp(ji,jj,jk) = tr(ji,jj,jk,jppic,Krhs) / rfact2/ (tr(ji,jj,jk,jppic,Kbb) + rtrn)
374        tr(ji,jj,jk,jpnpi,Krhs) = tr(ji,jj,jk,jpnpi,Krhs) + zproptot * texcretp
375        tr(ji,jj,jk,jpppi,Krhs) = tr(ji,jj,jk,jpppi,Krhs) + zpropo4p(ji,jj,jk) * texcretp   &
376        &                     + zprodopp(ji,jj,jk) * texcretp
377        tr(ji,jj,jk,jppfe,Krhs) = tr(ji,jj,jk,jppfe,Krhs) + zprofep(ji,jj,jk) * texcretp
378        tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) + zprorcad(ji,jj,jk) * texcretd   &
379           &                  - zpsino3 * zpronewd(ji,jj,jk) - zpsinh4 * zproregd(ji,jj,jk)   &
380           &                  - zrespd(ji,jj,jk) 
381        zcroissd(ji,jj,jk) = tr(ji,jj,jk,jpdia,Krhs) / rfact2 / (tr(ji,jj,jk,jpdia,Kbb) + rtrn)
382        tr(ji,jj,jk,jpndi,Krhs) = tr(ji,jj,jk,jpndi,Krhs) + zprodtot * texcretd
383        tr(ji,jj,jk,jppdi,Krhs) = tr(ji,jj,jk,jppdi,Krhs) + zpropo4d(ji,jj,jk) * texcretd   &
384        &                     + zprodopd(ji,jj,jk) * texcretd
385        tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) + zprofed(ji,jj,jk) * texcretd
386        tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd
387        tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
388        &                     + excretp * zprorcap(ji,jj,jk)
389        tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + excretd * zprodtot + excretn * zprontot   &
390        &                     + excretp * zproptot
391        tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + excretd * zpropo4d(ji,jj,jk) + excretn * zpropo4n(ji,jj,jk)   &
392        &    - texcretn * zprodopn(ji,jj,jk) - texcretd * zprodopd(ji,jj,jk) + excretp * zpropo4p(ji,jj,jk)     &
393        &    - texcretp * zprodopp(ji,jj,jk)
394        tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + o2ut * ( zproregn(ji,jj,jk) + zproregd(ji,jj,jk)   &
395           &                + zproregp(ji,jj,jk) ) + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk)           &
396           &                + zpronewd(ji,jj,jk) + zpronewp(ji,jj,jk) )   &
397           &                - o2ut * ( zrespn(ji,jj,jk) + zrespp(ji,jj,jk) + zrespd(ji,jj,jk) )
398        zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
399        tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zfeup
400        tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
401        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) - zprorcap(ji,jj,jk)  &
402        &                     + zpsino3 * zpronewn(ji,jj,jk) + zpsinh4 * zproregn(ji,jj,jk)   &
403        &                     + zpsino3 * zpronewp(ji,jj,jk) + zpsinh4 * zproregp(ji,jj,jk)   &
404        &                     + zpsino3 * zpronewd(ji,jj,jk) + zpsinh4 * zproregd(ji,jj,jk)  &
405        &                     + zrespn(ji,jj,jk) + zrespd(ji,jj,jk) + zrespp(ji,jj,jk) 
406        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk)  &
407        &                     + zpronewp(ji,jj,jk) ) - rno3 * ( zproregn(ji,jj,jk) + zproregd(ji,jj,jk)     &
408        &                     + zproregp(ji,jj,jk) ) 
409      END_3D
410     !
411     IF( ln_ligand ) THEN
412         zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp             
413         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
414           zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) + excretp * zprorcap(ji,jj,jk)
415           zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
416           tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet
417           zpligprod1(ji,jj,jk) = zdocprod * ldocp
418           zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet
419         END_3D
420     ENDIF
421
422
423     ! Total primary production per year
424
425    ! Total primary production per year
426    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
427      & tpp = glob_sum( 'p5zprod', ( zprorcan(:,:,:) + zprorcad(:,:,:) + zprorcap(:,:,:) ) * cvol(:,:,:) )
428
429    IF( lk_iomput .AND.  knt == nrdttrc ) THEN
430       zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
431       !
432       CALL iom_put( "PPPHYP"  , zprorcap(:,:,:) * zfact * tmask(:,:,:)   ) ! primary production by picophyto
433       CALL iom_put( "PPPHYN"  , zprorcan(:,:,:) * zfact * tmask(:,:,:) )  ! primary production by nanophyto
434       CALL iom_put( "PPPHYD"  , zprorcad(:,:,:) * zfact * tmask(:,:,:)   ) ! primary production by diatomes
435       CALL iom_put( "PPNEWN"  , zpronewp(:,:,:) * zfact * tmask(:,:,:)    ) ! new primary production by picophyto
436       CALL iom_put( "PPNEWN"  , zpronewn(:,:,:) * zfact * tmask(:,:,:)    ) ! new primary production by nanophyto
437       CALL iom_put( "PPNEWD"  , zpronewd(:,:,:) * zfact * tmask(:,:,:)   ) ! new primary production by diatomes
438       CALL iom_put( "PBSi"    , zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)  ) ! biogenic silica production
439       CALL iom_put( "PFeP"    , zprofep(:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by picophyto
440       CALL iom_put( "PFeN"    , zprofen(:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by nanophyto
441       CALL iom_put( "PFeD"    , zprofed(:,:,:) * zfact * tmask(:,:,:)  ) ! biogenic iron production by  diatomes
442       IF( ln_ligand ) THEN
443         CALL iom_put( "LPRODP"  , zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:) )
444         CALL iom_put( "LDETP"   , zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:) )
445       ENDIF
446       CALL iom_put( "Mumax"   , zprmaxn(:,:,:) * tmask(:,:,:)  ) ! Maximum growth rate
447       CALL iom_put( "MuP"     , zprpic(:,:,:) * xlimpic(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for picophyto
448       CALL iom_put( "MuN"     , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for nanophyto
449       CALL iom_put( "MuD"     , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for diatoms
450       CALL iom_put( "LPlight" , zprpic(:,:,:) / (zprmaxp(:,:,:) + rtrn) * tmask(:,:,:)  )  ! light limitation term
451       CALL iom_put( "LNlight" , zprbio(:,:,:) / (zprmaxn(:,:,:) + rtrn) * tmask(:,:,:)  )  ! light limitation term
452       CALL iom_put( "LDlight" , zprdia(:,:,:) / (zprmaxd(:,:,:) + rtrn) * tmask(:,:,:)   )
453       CALL iom_put( "MunetP"  , zcroissp(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for picophyto
454       CALL iom_put( "MunetN"  , zcroissn(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for nanophyto
455       CALL iom_put( "MunetD"  , zcroissd(:,:,:) * tmask(:,:,:) ) ! Realized growth rate for diatoms
456       CALL iom_put( "TPP"     , ( zprorcap(:,:,:) + zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  )  ! total primary production
457       CALL iom_put( "TPNEW"   , ( zpronewp(:,:,:) + zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ) ! total new production
458       CALL iom_put( "TPBFE"   , ( zprofep (:,:,:) + zprofen (:,:,:) + zprofed (:,:,:) ) * zfact * tmask(:,:,:)  )  ! total biogenic iron production
459       CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
460     ENDIF
461
462      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
463         WRITE(charout, FMT="('prod')")
464         CALL prt_ctl_info( charout, cdcomp = 'top' )
465         CALL prt_ctl(tab4d_1=CASTWP(tr(:,:,:,:,Krhs)), mask1=tmask, clinfo=ctrcnm)
466      ENDIF
467      !
468      IF( ln_timing )   CALL timing_stop('p5z_prod')
469      !
470   END SUBROUTINE p5z_prod
471
472
473   SUBROUTINE p5z_prod_init
474      !!----------------------------------------------------------------------
475      !!                  ***  ROUTINE p5z_prod_init  ***
476      !!
477      !! ** Purpose :   Initialization of phytoplankton production parameters
478      !!
479      !! ** Method  :   Read the nampisprod namelist and check the parameters
480      !!      called at the first timestep (nittrc000)
481      !!
482      !! ** input   :   Namelist nampisprod
483      !!----------------------------------------------------------------------
484      INTEGER :: ios                 ! Local integer output status for namelist read
485      !!
486      NAMELIST/namp5zprod/ pislopen, pislopep, pisloped, excretn, excretp, excretd,     &
487         &                 thetannm, thetanpm, thetandm, chlcmin, grosip, bresp, xadap
488      !!----------------------------------------------------------------------
489
490      READ  ( numnatp_ref, namp5zprod, IOSTAT = ios, ERR = 901)
491901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zprod in reference namelist' )
492
493      READ  ( numnatp_cfg, namp5zprod, IOSTAT = ios, ERR = 902 )
494902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zprod in configuration namelist' )
495      IF(lwm) WRITE ( numonp, namp5zprod )
496
497      IF(lwp) THEN                         ! control print
498         WRITE(numout,*) ' '
499         WRITE(numout,*) ' Namelist parameters for phytoplankton growth, namp5zprod'
500         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
501         WRITE(numout,*) '    mean Si/C ratio                           grosip       =', grosip
502         WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
503         WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
504         WRITE(numout,*) '    P-I slope  for picophytoplankton          pislopep     =', pislopep
505         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
506         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
507         WRITE(numout,*) '    excretion ratio of picophytoplankton      excretp      =', excretp
508         WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd
509         WRITE(numout,*) '    basal respiration in phytoplankton        bresp        =', bresp
510         WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
511         WRITE(numout,*) '    Minimum Chl/N in nanophytoplankton        thetannm     =', thetannm
512         WRITE(numout,*) '    Minimum Chl/N in picophytoplankton        thetanpm     =', thetanpm
513         WRITE(numout,*) '    Minimum Chl/N in diatoms                  thetandm     =', thetandm
514      ENDIF
515      !
516      r1_rday   = 1._wp / rday 
517      texcretn  = 1._wp - excretn
518      texcretp  = 1._wp - excretp
519      texcretd  = 1._wp - excretd
520      tpp       = 0._wp
521      !
522   END SUBROUTINE p5z_prod_init
523
524
525   INTEGER FUNCTION p5z_prod_alloc()
526      !!----------------------------------------------------------------------
527      !!                     ***  ROUTINE p5z_prod_alloc  ***
528      !!----------------------------------------------------------------------
529      ALLOCATE( zdaylen(jpi,jpj), STAT = p5z_prod_alloc )
530      !
531      IF( p5z_prod_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p5z_prod_alloc : failed to allocate arrays.' )
532      !
533   END FUNCTION p5z_prod_alloc
534   !!======================================================================
535END MODULE p5zprod
Note: See TracBrowser for help on using the repository browser.