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

Last change on this file since 10973 was 10973, checked in by cetlod, 19 months ago

trunk : bugfix in Lability of big particles, see ticket #2278

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