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/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 11317

Last change on this file since 11317 was 11317, checked in by smasson, 5 years ago

dev_r10984_HPC-13 : improve error handling, see #2307 and #2285

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