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.
p4zpoc.F90 in NEMO/trunk/src/TOP/PISCES/P4Z – NEMO

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

Last change on this file since 10069 was 10069, checked in by nicolasmartin, 6 years ago

Fix mistakes of previous commit on SVN keywords property

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