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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 11463

Last change on this file since 11463 was 10975, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Finish converting all TOP routines and knock-on effects of these conversions. Fully SETTE tested (SETTE tests 1-6 and 9). This completes the first stage conversion of TRA and TOP but need to revisit and pass ts and tr arrays through the argument lists where appropriate.

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