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/2020/dev_r12563_ASINTER-06_ABL_improvement/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 13891

Last change on this file since 13891 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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