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.
p5zpoc.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z – NEMO

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90 @ 6966

Last change on this file since 6966 was 6966, checked in by aumont, 7 years ago

various bug fixes and updates on carbon chemistry

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