source: NEMO/trunk/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 32.5 KB
Line 
1MODULE p4zpoc
2   !!======================================================================
3   !!                         ***  MODULE p4zpoc  ***
4   !! TOP :   PISCES Compute remineralization of organic particles
5   !!=========================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Quota model for iron
9   !!             3.6  !  2016-03  (O. Aumont) Quota model and diverse
10   !!----------------------------------------------------------------------
11   !!   p4z_poc       :  Compute remineralization/dissolution of organic compounds
12   !!   p4z_poc_init  :  Initialisation of parameters for remineralisation
13   !!----------------------------------------------------------------------
14   USE oce_trc         !  shared variables between ocean and passive tracers
15   USE trc             !  passive tracers common variables
16   USE sms_pisces      !  PISCES Source Minus Sink variables
17   USE prtctl_trc      !  print control for debugging
18   USE iom             !  I/O manager
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   p4z_poc         ! called in p4zbio.F90
24   PUBLIC   p4z_poc_init    ! called in trcsms_pisces.F90
25   PUBLIC   alngam          !
26   PUBLIC   gamain          !
27
28   REAL(wp), PUBLIC ::   xremip     !: remineralisation rate of DOC
29   REAL(wp), PUBLIC ::   xremipc    !: remineralisation rate of DOC
30   REAL(wp), PUBLIC ::   xremipn    !: remineralisation rate of DON
31   REAL(wp), PUBLIC ::   xremipp    !: remineralisation rate of DOP
32   INTEGER , PUBLIC ::   jcpoc      !: number of lability classes
33   REAL(wp), PUBLIC ::   rshape     !: shape factor of the gamma distribution
34
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)       ::   alphan, reminp   !:
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   alphap           !:
37
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
44   !! $Id$
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE p4z_poc( kt, knt, Kbb, Kmm, Krhs )
50      !!---------------------------------------------------------------------
51      !!                     ***  ROUTINE p4z_poc  ***
52      !!
53      !! ** Purpose :   Compute remineralization of organic particles
54      !!
55      !! ** Method  : - ???
56      !!---------------------------------------------------------------------
57      INTEGER, INTENT(in) ::   kt, knt         ! ocean time step and ???
58      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
59      !
60      INTEGER  ::   ji, jj, jk, jn
61      REAL(wp) ::   zremip, zremig, zdep, zorem, zorem2, zofer
62      REAL(wp) ::   zopon, zopop, zopoc, zopoc2, zopon2, zopop2
63      REAL(wp) ::   zsizek, zsizek1, alphat, remint, solgoc, zpoc
64      REAL(wp) ::   zofer2, zofer3
65      REAL(wp) ::   zrfact2
66      CHARACTER (len=25) :: charout
67      REAL(wp), DIMENSION(jpi,jpj  )   :: totprod, totthick, totcons 
68      REAL(wp), DIMENSION(jpi,jpj,jpk)   :: zremipoc, zremigoc, zorem3, ztremint, zfolimi
69      REAL(wp), DIMENSION(jpi,jpj,jpk,jcpoc) :: alphag
70      !!---------------------------------------------------------------------
71      !
72      IF( ln_timing )  CALL timing_start('p4z_poc')
73      !
74      ! Initialization of local variables
75      ! ---------------------------------
76
77      ! Here we compute the GOC -> POC rate due to the shrinking
78      ! of the fecal pellets/aggregates as a result of bacterial
79      ! solubilization
80      ! This is based on a fractal dimension of 2.56 and a spectral
81      ! slope of -3.6 (identical to what is used in p4zsink to compute
82      ! aggregation
83      solgoc = 0.04/ 2.56 * 1./ ( 1.-50**(-0.04) )
84
85      ! Initialisation of temprary arrys
86      IF( ln_p4z ) THEN
87         zremipoc(:,:,:) = xremip
88         zremigoc(:,:,:) = xremip
89      ELSE    ! ln_p5z
90         zremipoc(:,:,:) = xremipc
91         zremigoc(:,:,:) = xremipc
92      ENDIF
93      zorem3(:,:,:)   = 0.
94      orem  (:,:,:)   = 0.
95      ztremint(:,:,:) = 0.
96      zfolimi (:,:,:) = 0.
97
98      DO jn = 1, jcpoc
99        alphag(:,:,:,jn) = alphan(jn)
100        alphap(:,:,:,jn) = alphan(jn)
101      END DO
102
103     ! -----------------------------------------------------------------------
104     ! Lability parameterization. This is the big particles part (GOC)
105     ! This lability parameterization can be activated only with the standard
106     ! particle scheme. Does not work with Kriest parameterization.
107     ! -----------------------------------------------------------------------
108     ztremint(:,:,:) = zremigoc(:,:,:)
109     DO_3D_11_11( 2, jpkm1 )
110        IF (tmask(ji,jj,jk) == 1.) THEN
111          zdep = hmld(ji,jj)
112          !
113          ! In the case of GOC, lability is constant in the mixed layer
114          ! It is computed only below the mixed layer depth
115          ! ------------------------------------------------------------
116          !
117          IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
118            alphat = 0.
119            remint = 0.
120            !
121            zsizek1  = e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
122            zsizek = e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
123            !
124            IF ( gdept(ji,jj,jk-1,Kmm) <= zdep ) THEN
125              !
126              ! The first level just below the mixed layer needs a
127              ! specific treatment because lability is supposed constant
128              ! everywhere within the mixed layer. This means that
129              ! change in lability in the bottom part of the previous cell
130              ! should not be computed
131              ! ----------------------------------------------------------
132              !
133              ! POC concentration is computed using the lagrangian
134              ! framework. It is only used for the lability param
135              zpoc = tr(ji,jj,jk-1,jpgoc,Kbb) + consgoc(ji,jj,jk) * rday / rfact2               &
136              &   * e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn)
137              zpoc = MAX(0., zpoc)
138              !
139              DO jn = 1, jcpoc
140                 !
141                 ! Lagrangian based algorithm. The fraction of each
142                 ! lability class is computed starting from the previous
143                 ! level
144                 ! -----------------------------------------------------
145                 !
146                 ! the concentration of each lability class is calculated
147                 ! as the sum of the different sources and sinks
148                 ! Please note that production of new GOC experiences
149                 ! degradation
150                 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
151                 &   + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn)             &
152                 &   * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2 
153                 alphat = alphat + alphag(ji,jj,jk,jn)
154                 remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
155              END DO
156            ELSE
157              !
158              ! standard algorithm in the rest of the water column
159              ! See the comments in the previous block.
160              ! ---------------------------------------------------
161              !
162              zpoc = tr(ji,jj,jk-1,jpgoc,Kbb) + consgoc(ji,jj,jk-1) * rday / rfact2               &
163              &   * e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
164              &   * rday / rfact2 * e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn)
165              zpoc = max(0., zpoc)
166              !
167              DO jn = 1, jcpoc
168                 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek              &
169                 &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
170                 &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
171                 &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn) 
172                 alphat = alphat + alphag(ji,jj,jk,jn)
173                 remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
174              END DO
175            ENDIF
176            !
177            DO jn = 1, jcpoc
178               ! The contribution of each lability class at the current
179               ! level is computed
180               alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn)
181            END DO
182            ! Computation of the mean remineralisation rate
183            ztremint(ji,jj,jk) =  MAX(0., remint / ( alphat + rtrn) )
184            !
185          ENDIF
186        ENDIF
187     END_3D
188
189      IF( ln_p4z ) THEN   ;   zremigoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
190      ELSE                ;   zremigoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
191      ENDIF
192
193      IF( ln_p4z ) THEN
194         DO_3D_11_11( 1, jpkm1 )
195            ! POC disaggregation by turbulence and bacterial activity.
196            ! --------------------------------------------------------
197            zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
198            zorem2  = zremig * tr(ji,jj,jk,jpgoc,Kbb)
199            orem(ji,jj,jk)      = zorem2
200            zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
201            zofer2 = zremig * tr(ji,jj,jk,jpbfe,Kbb)
202            zofer3 = zremig * solgoc * tr(ji,jj,jk,jpbfe,Kbb)
203
204            ! -------------------------------------
205            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
206            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zorem2 - zorem3(ji,jj,jk)
207            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zofer3
208            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 - zofer3
209            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem2
210            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
211            zfolimi(ji,jj,jk)   = zofer2
212         END_3D
213      ELSE
214         DO_3D_11_11( 1, jpkm1 )
215             ! POC disaggregation by turbulence and bacterial activity.
216            ! --------------------------------------------------------
217            zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
218            zopoc2 = zremig  * tr(ji,jj,jk,jpgoc,Kbb)
219            orem(ji,jj,jk) = zopoc2
220            zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
221            zopon2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpgon,Kbb)
222            zopop2 = xremipp / xremipc * zremig * tr(ji,jj,jk,jpgop,Kbb)
223            zofer2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpbfe,Kbb)
224
225            ! -------------------------------------
226            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
227            tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + solgoc * zopon2 
228            tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + solgoc * zopop2
229            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + solgoc * zofer2
230            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc2
231            tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon2
232            tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop2
233            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
234            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zopoc2 - zorem3(ji,jj,jk)
235            tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) - zopon2 * (1. + solgoc)
236            tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) - zopop2 * (1. + solgoc)
237            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 * (1. + solgoc)
238            zfolimi(ji,jj,jk)   = zofer2
239         END_3D
240      ENDIF
241
242     IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
243        WRITE(charout, FMT="('poc1')")
244        CALL prt_ctl_trc_info(charout)
245        CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
246     ENDIF
247
248     ! ------------------------------------------------------------------
249     ! Lability parameterization for the small OM particles. This param
250     ! is based on the same theoretical background as the big particles.
251     ! However, because of its low sinking speed, lability is not supposed
252     ! to be equal to its initial value (the value of the freshly produced
253     ! organic matter). It is however uniform in the mixed layer.
254     ! -------------------------------------------------------------------
255     !
256     totprod (:,:) = 0.
257     totthick(:,:) = 0.
258     totcons (:,:) = 0.
259     ! intregrated production and consumption of POC in the mixed layer
260     ! ----------------------------------------------------------------
261     !
262     DO_3D_11_11( 1, jpkm1 )
263        zdep = hmld(ji,jj)
264        IF (tmask(ji,jj,jk) == 1. .AND. gdept(ji,jj,jk,Kmm) <= zdep ) THEN
265          totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2
266          ! The temperature effect is included here
267          totthick(ji,jj) = totthick(ji,jj) + e3t(ji,jj,jk,Kmm)* tgfunc(ji,jj,jk)
268          totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2    &
269          &                / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn )
270        ENDIF
271     END_3D
272
273     ! Computation of the lability spectrum in the mixed layer. In the mixed
274     ! layer, this spectrum is supposed to be uniform.
275     ! ---------------------------------------------------------------------
276     ztremint(:,:,:) = zremipoc(:,:,:)
277     DO_3D_11_11( 1, jpkm1 )
278        IF (tmask(ji,jj,jk) == 1.) THEN
279          zdep = hmld(ji,jj)
280          alphat = 0.0
281          remint = 0.0
282          IF( gdept(ji,jj,jk,Kmm) <= zdep ) THEN
283             DO jn = 1, jcpoc
284                ! For each lability class, the system is supposed to be
285                ! at equilibrium: Prod - Sink - w alphap = 0.
286                alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn)    &
287                &                     * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
288                alphat = alphat + alphap(ji,jj,jk,jn)
289             END DO
290             DO jn = 1, jcpoc
291                alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
292                remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
293             END DO
294             ! Mean remineralization rate in the mixed layer
295             ztremint(ji,jj,jk) =  MAX( 0., remint )
296          ENDIF
297        ENDIF
298     END_3D
299     !
300     IF( ln_p4z ) THEN   ;  zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
301     ELSE                ;  zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
302     ENDIF
303
304     ! -----------------------------------------------------------------------
305     ! The lability parameterization is used here. The code is here
306     ! almost identical to what is done for big particles. The only difference
307     ! is that an additional source from GOC to POC is included. This means
308     ! that since we need the lability spectrum of GOC, GOC spectrum
309     ! should be determined before.
310     ! -----------------------------------------------------------------------
311     !
312     DO_3D_11_11( 2, jpkm1 )
313        IF (tmask(ji,jj,jk) == 1.) THEN
314          zdep = hmld(ji,jj)
315          IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
316            alphat = 0.
317            remint = 0.
318            !
319            ! the scale factors are corrected with temperature
320            zsizek1  = e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
321            zsizek = e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
322            !
323            ! Special treatment of the level just below the MXL
324            ! See the comments in the GOC section
325            ! ---------------------------------------------------
326            !
327            IF ( gdept(ji,jj,jk-1,Kmm) <= zdep ) THEN
328              !
329              ! Computation of the POC concentration using the
330              ! lagrangian algorithm
331              zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk) * rday / rfact2               &
332              &   * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
333              zpoc = max(0., zpoc)
334              !
335              DO jn = 1, jcpoc
336                 ! computation of the lability spectrum applying the
337                 ! different sources and sinks
338                 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc  &
339                 &   + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
340                 &   / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn)     &
341                 &   * zsizek ) )
342                 alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
343                 alphat = alphat + alphap(ji,jj,jk,jn)
344              END DO
345            ELSE
346              !
347              ! Lability parameterization for the interior of the ocean
348              ! This is very similar to what is done in the previous
349              ! block
350              ! --------------------------------------------------------
351              !
352              zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk-1) * rday / rfact2               &
353              &   * e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
354              &   * rday / rfact2 * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
355              zpoc = max(0., zpoc)
356              !
357              DO jn = 1, jcpoc
358                 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)                       &
359                 &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn)             & 
360                 &   + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn)      &
361                 &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn)  &
362                 &   * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk)                 &
363                 &   * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1.   &
364                 &   - exp( -reminp(jn) * zsizek ) )
365                 alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
366                 alphat = alphat + alphap(ji,jj,jk,jn)
367              END DO
368            ENDIF
369            ! Normalization of the lability spectrum so that the
370            ! integral is equal to 1
371            DO jn = 1, jcpoc
372               alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
373               remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
374            END DO
375            ! Mean remineralization rate in the water column
376            ztremint(ji,jj,jk) =  MAX( 0., remint )
377          ENDIF
378        ENDIF
379     END_3D
380
381     IF( ln_p4z ) THEN   ;   zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
382     ELSE                ;   zremipoc(:,:,:) = MIN( xremipc, ztremint(:,:,:) )
383     ENDIF
384
385     IF( ln_p4z ) THEN
386         DO_3D_11_11( 1, jpkm1 )
387            IF (tmask(ji,jj,jk) == 1.) THEN
388              ! POC disaggregation by turbulence and bacterial activity.
389              ! --------------------------------------------------------
390              zremip          = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
391              zorem           = zremip * tr(ji,jj,jk,jppoc,Kbb)
392              zofer           = zremip * tr(ji,jj,jk,jpsfe,Kbb)
393
394              tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem
395              orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
396              tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer
397              tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zorem
398              tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
399              zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
400            ENDIF
401         END_3D
402     ELSE
403       DO_3D_11_11( 1, jpkm1 )
404          ! POC disaggregation by turbulence and bacterial activity.
405          ! --------------------------------------------------------
406          zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
407          zopoc  = zremip * tr(ji,jj,jk,jppoc,Kbb)
408          orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
409          zopon  = xremipn / xremipc * zremip * tr(ji,jj,jk,jppon,Kbb)
410          zopop  = xremipp / xremipc * zremip * tr(ji,jj,jk,jppop,Kbb)
411          zofer  = xremipn / xremipc * zremip * tr(ji,jj,jk,jpsfe,Kbb)
412
413          tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zopoc
414          tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zopon
415          tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) - zopop
416          tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
417          tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc
418          tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon 
419          tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop 
420          tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer 
421          zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
422       END_3D
423     ENDIF
424
425     IF( lk_iomput ) THEN
426        IF( knt == nrdttrc ) THEN
427          zrfact2 = 1.e3 * rfact2r
428          CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
429          CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
430          CALL iom_put( "REMINF" , zfolimi(:,:,:)  * tmask(:,:,:)  * 1.e+9 * zrfact2 )  ! Remineralisation rate
431        ENDIF
432     ENDIF
433
434      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
435         WRITE(charout, FMT="('poc2')")
436         CALL prt_ctl_trc_info(charout)
437         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
438      ENDIF
439      !
440      !
441      IF( ln_timing )   CALL timing_stop('p4z_poc')
442      !
443   END SUBROUTINE p4z_poc
444
445
446   SUBROUTINE p4z_poc_init
447      !!----------------------------------------------------------------------
448      !!                  ***  ROUTINE p4z_poc_init  ***
449      !!
450      !! ** Purpose :   Initialization of remineralization parameters
451      !!
452      !! ** Method  :   Read the nampispoc namelist and check the parameters
453      !!              called at the first timestep
454      !!
455      !! ** input   :   Namelist nampispoc
456      !!----------------------------------------------------------------------
457      INTEGER ::   jn            ! dummy loop index
458      INTEGER ::   ios, ifault   ! Local integer
459      REAL(wp)::   remindelta, reminup, remindown
460      !!
461      NAMELIST/nampispoc/ xremip , jcpoc  , rshape,  &
462         &                xremipc, xremipn, xremipp
463      !!----------------------------------------------------------------------
464      !
465      IF(lwp) THEN
466         WRITE(numout,*)
467         WRITE(numout,*) 'p4z_poc_init : Initialization of remineralization parameters'
468         WRITE(numout,*) '~~~~~~~~~~~~'
469      ENDIF
470      !
471      READ  ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
472901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampispoc in reference namelist' )
473      READ  ( numnatp_cfg, nampispoc, IOSTAT = ios, ERR = 902 )
474902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampispoc in configuration namelist' )
475      IF(lwm) WRITE( numonp, nampispoc )
476
477      IF(lwp) THEN                         ! control print
478         WRITE(numout,*) '   Namelist : nampispoc'
479         IF( ln_p4z ) THEN
480            WRITE(numout,*) '      remineralisation rate of POC              xremip    =', xremip
481         ELSE
482            WRITE(numout,*) '      remineralisation rate of POC              xremipc   =', xremipc
483            WRITE(numout,*) '      remineralisation rate of PON              xremipn   =', xremipn
484            WRITE(numout,*) '      remineralisation rate of POP              xremipp   =', xremipp
485         ENDIF
486         WRITE(numout,*) '      Number of lability classes for POC        jcpoc     =', jcpoc
487         WRITE(numout,*) '      Shape factor of the gamma distribution    rshape    =', rshape
488      ENDIF
489      !
490      ! Discretization along the lability space
491      ! ---------------------------------------
492      !
493      ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
494      !
495      IF (jcpoc > 1) THEN
496         !
497         remindelta = LOG(4. * 1000. ) / REAL(jcpoc-1, wp)
498         reminup = 1./ 400. * EXP(remindelta)
499         !
500         ! Discretization based on incomplete gamma functions
501         ! As incomplete gamma functions are not available in standard
502         ! fortran 95, they have been coded as functions in this module (gamain)
503         ! ---------------------------------------------------------------------
504         !
505         alphan(1) = gamain(reminup, rshape, ifault)
506         reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1)
507         DO jn = 2, jcpoc-1
508            reminup = 1./ 400. * EXP( REAL(jn, wp) * remindelta)
509            remindown = 1. / 400. * EXP( REAL(jn-1, wp) * remindelta)
510            alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault)
511            reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault)
512            reminp(jn) = reminp(jn) * xremip / alphan(jn)
513         END DO
514         remindown = 1. / 400. * EXP( REAL(jcpoc-1, wp) * remindelta)
515         alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault)
516         reminp(jcpoc) = 1.0 - gamain(remindown, rshape+1.0, ifault)
517         reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
518
519      ELSE
520         alphan(jcpoc) = 1.
521         reminp(jcpoc) = xremip
522      ENDIF
523
524      DO jn = 1, jcpoc
525         alphap(:,:,:,jn) = alphan(jn)
526      END DO
527
528   END SUBROUTINE p4z_poc_init
529
530
531   REAL FUNCTION alngam( xvalue, ifault )
532      !*****************************************************************************80
533      !
534      !! ALNGAM computes the logarithm of the gamma function.
535      !
536      !  Modified:    13 January 2008
537      !
538      !  Author  :    Allan Macleod
539      !               FORTRAN90 version by John Burkardt
540      !
541      !  Reference:
542      !    Allan Macleod, Algorithm AS 245,
543      !    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
544      !    Applied Statistics,
545      !    Volume 38, Number 2, 1989, pages 397-402.
546      !
547      !  Parameters:
548      !
549      !    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
550      !
551      !    Output, integer ( kind = 4 ) IFAULT, error flag.
552      !    0, no error occurred.
553      !    1, XVALUE is less than or equal to 0.
554      !    2, XVALUE is too big.
555      !
556      !    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
557      !*****************************************************************************80
558  implicit none
559
560  real(wp), parameter :: alr2pi = 0.918938533204673E+00
561  integer:: ifault
562  real(wp), dimension ( 9 ) :: r1 = (/ &
563    -2.66685511495E+00, &
564    -24.4387534237E+00, &
565    -21.9698958928E+00, &
566     11.1667541262E+00, &
567     3.13060547623E+00, &
568     0.607771387771E+00, &
569     11.9400905721E+00, &
570     31.4690115749E+00, &
571     15.2346874070E+00 /)
572  real(wp), dimension ( 9 ) :: r2 = (/ &
573    -78.3359299449E+00, &
574    -142.046296688E+00, &
575     137.519416416E+00, &
576     78.6994924154E+00, &
577     4.16438922228E+00, &
578     47.0668766060E+00, &
579     313.399215894E+00, &
580     263.505074721E+00, &
581     43.3400022514E+00 /)
582  real(wp), dimension ( 9 ) :: r3 = (/ &
583    -2.12159572323E+05, &
584     2.30661510616E+05, &
585     2.74647644705E+04, &
586    -4.02621119975E+04, &
587    -2.29660729780E+03, &
588    -1.16328495004E+05, &
589    -1.46025937511E+05, &
590    -2.42357409629E+04, &
591    -5.70691009324E+02 /)
592  real(wp), dimension ( 5 ) :: r4 = (/ &
593     0.279195317918525E+00, &
594     0.4917317610505968E+00, &
595     0.0692910599291889E+00, &
596     3.350343815022304E+00, &
597     6.012459259764103E+00 /)
598  real (wp) :: x
599  real (wp) :: x1
600  real (wp) :: x2
601  real (wp), parameter :: xlge = 5.10E+05
602  real (wp), parameter :: xlgst = 1.0E+30
603  real (wp) :: xvalue
604  real (wp) :: y
605
606  x = xvalue
607  alngam = 0.0E+00
608!
609!  Check the input.
610!
611  if ( xlgst <= x ) then
612    ifault = 2
613    return
614  end if
615  if ( x <= 0.0E+00 ) then
616    ifault = 1
617    return
618  end if
619
620  ifault = 0
621!
622!  Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
623!
624  if ( x < 1.5E+00 ) then
625
626    if ( x < 0.5E+00 ) then
627      alngam = - log ( x )
628      y = x + 1.0E+00
629!
630!  Test whether X < machine epsilon.
631!
632      if ( y == 1.0E+00 ) then
633        return
634      end if
635
636    else
637
638      alngam = 0.0E+00
639      y = x
640      x = ( x - 0.5E+00 ) - 0.5E+00
641
642    end if
643
644    alngam = alngam + x * (((( &
645        r1(5)   * y &
646      + r1(4) ) * y &
647      + r1(3) ) * y &
648      + r1(2) ) * y &
649      + r1(1) ) / (((( &
650                  y &
651      + r1(9) ) * y &
652      + r1(8) ) * y &
653      + r1(7) ) * y &
654      + r1(6) )
655
656    return
657
658  end if
659!
660!  Calculation for 1.5 <= X < 4.0.
661!
662  if ( x < 4.0E+00 ) then
663
664    y = ( x - 1.0E+00 ) - 1.0E+00
665
666    alngam = y * (((( &
667        r2(5)   * x &
668      + r2(4) ) * x &
669      + r2(3) ) * x &
670      + r2(2) ) * x &
671      + r2(1) ) / (((( &
672                  x &
673      + r2(9) ) * x &
674      + r2(8) ) * x &
675      + r2(7) ) * x &
676      + r2(6) )
677!
678!  Calculation for 4.0 <= X < 12.0.
679!
680  else if ( x < 12.0E+00 ) then
681
682    alngam = (((( &
683        r3(5)   * x &
684      + r3(4) ) * x &
685      + r3(3) ) * x &
686      + r3(2) ) * x &
687      + r3(1) ) / (((( &
688                  x &
689      + r3(9) ) * x &
690      + r3(8) ) * x &
691      + r3(7) ) * x &
692      + r3(6) )
693!
694!  Calculation for 12.0 <= X.
695!
696  else
697
698    y = log ( x )
699    alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
700
701    if ( x <= xlge ) then
702
703      x1 = 1.0E+00 / x
704      x2 = x1 * x1
705
706      alngam = alngam + x1 * ( ( &
707             r4(3)   * &
708        x2 + r4(2) ) * &
709        x2 + r4(1) ) / ( ( &
710        x2 + r4(5) ) * &
711        x2 + r4(4) )
712
713    end if
714
715   end if
716
717   END FUNCTION alngam
718
719
720   REAL FUNCTION gamain( x, p, ifault )
721!*****************************************************************************80
722!
723!! GAMAIN computes the incomplete gamma ratio.
724!
725!  Discussion:
726!
727!    A series expansion is used if P > X or X <= 1.  Otherwise, a
728!    continued fraction approximation is used.
729!
730!  Modified:
731!
732!    17 January 2008
733!
734!  Author:
735!
736!    G Bhattacharjee
737!    FORTRAN90 version by John Burkardt
738!
739!  Reference:
740!
741!    G Bhattacharjee,
742!    Algorithm AS 32:
743!    The Incomplete Gamma Integral,
744!    Applied Statistics,
745!    Volume 19, Number 3, 1970, pages 285-287.
746!
747!  Parameters:
748!
749!    Input, real ( kind = 8 ) X, P, the parameters of the incomplete
750!    gamma ratio.  0 <= X, and 0 < P.
751!
752!    Output, integer ( kind = 4 ) IFAULT, error flag.
753!    0, no errors.
754!    1, P <= 0.
755!    2, X < 0.
756!    3, underflow.
757!    4, error return from the Log Gamma routine.
758!
759!    Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
760!    gamma ratio.
761!
762  implicit none
763
764  real (wp) a
765  real (wp), parameter :: acu = 1.0E-08
766  real (wp) an
767  real (wp) arg
768  real (wp) b
769  real (wp) dif
770  real (wp) factor
771  real (wp) g
772  real (wp) gin
773  integer i
774  integer ifault
775  real (wp), parameter :: oflo = 1.0E+37
776  real (wp) p
777  real (wp) pn(6)
778  real (wp) rn
779  real (wp) term
780  real (wp), parameter :: uflo = 1.0E-37
781  real (wp) x
782!
783!  Check the input.
784!
785  if ( p <= 0.0E+00 ) then
786    ifault = 1
787    gamain = 0.0E+00
788    return
789  end if
790
791  if ( x < 0.0E+00 ) then
792    ifault = 2
793    gamain = 0.0E+00
794    return
795  end if
796
797  if ( x == 0.0E+00 ) then
798    ifault = 0
799    gamain = 0.0E+00
800    return
801  end if
802
803  g = alngam ( p, ifault )
804
805  if ( ifault /= 0 ) then
806    ifault = 4
807    gamain = 0.0E+00
808    return
809  end if
810
811  arg = p * log ( x ) - x - g
812
813  if ( arg < log ( uflo ) ) then
814    ifault = 3
815    gamain = 0.0E+00
816    return
817  end if
818
819  ifault = 0
820  factor = exp ( arg )
821!
822!  Calculation by series expansion.
823!
824  if ( x <= 1.0E+00 .or. x < p ) then
825
826    gin = 1.0E+00
827    term = 1.0E+00
828    rn = p
829
830    do
831
832      rn = rn + 1.0E+00
833      term = term * x / rn
834      gin = gin + term
835
836      if ( term <= acu ) then
837        exit
838      end if
839
840    end do
841
842    gamain = gin * factor / p
843    return
844
845  end if
846!
847!  Calculation by continued fraction.
848!
849  a = 1.0E+00 - p
850  b = a + x + 1.0E+00
851  term = 0.0E+00
852
853  pn(1) = 1.0E+00
854  pn(2) = x
855  pn(3) = x + 1.0E+00
856  pn(4) = x * b
857
858  gin = pn(3) / pn(4)
859
860  do
861
862    a = a + 1.0E+00
863    b = b + 2.0E+00
864    term = term + 1.0E+00
865    an = a * term
866    do i = 1, 2
867      pn(i+4) = b * pn(i+2) - an * pn(i)
868    end do
869
870    if ( pn(6) /= 0.0E+00 ) then
871
872      rn = pn(5) / pn(6)
873      dif = abs ( gin - rn )
874!
875!  Absolute error tolerance satisfied?
876!
877      if ( dif <= acu ) then
878!
879!  Relative error tolerance satisfied?
880!
881        if ( dif <= acu * rn ) then
882          gamain = 1.0E+00 - factor * gin
883          exit
884        end if
885
886      end if
887
888      gin = rn
889
890    end if
891
892    do i = 1, 4
893      pn(i) = pn(i+2)
894    end do
895    if ( oflo <= abs ( pn(5) ) ) then
896
897      do i = 1, 4
898        pn(i) = pn(i) / oflo
899      end do
900
901    end if
902
903  end do
904
905  END FUNCTION gamain
906
907   !!======================================================================
908END MODULE p4zpoc
Note: See TracBrowser for help on using the repository browser.