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 @ 7190

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

various important bug fixes

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