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_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/PISCES/P4Z/p4zpoc.F90 @ 11671

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

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.

  • Property svn:keywords set to Id
File size: 34.0 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      READ  ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
500901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampispoc in reference namelist' )
501      READ  ( numnatp_cfg, nampispoc, IOSTAT = ios, ERR = 902 )
502902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampispoc in configuration namelist' )
503      IF(lwm) WRITE( numonp, nampispoc )
504
505      IF(lwp) THEN                         ! control print
506         WRITE(numout,*) '   Namelist : nampispoc'
507         IF( ln_p4z ) THEN
508            WRITE(numout,*) '      remineralisation rate of POC              xremip    =', xremip
509         ELSE
510            WRITE(numout,*) '      remineralisation rate of POC              xremipc   =', xremipc
511            WRITE(numout,*) '      remineralisation rate of PON              xremipn   =', xremipn
512            WRITE(numout,*) '      remineralisation rate of POP              xremipp   =', xremipp
513         ENDIF
514         WRITE(numout,*) '      Number of lability classes for POC        jcpoc     =', jcpoc
515         WRITE(numout,*) '      Shape factor of the gamma distribution    rshape    =', rshape
516      ENDIF
517      !
518      ! Discretization along the lability space
519      ! ---------------------------------------
520      !
521      ALLOCATE( alphan(jcpoc) , reminp(jcpoc) , alphap(jpi,jpj,jpk,jcpoc) )
522      !
523      IF (jcpoc > 1) THEN
524         !
525         remindelta = LOG(4. * 1000. ) / REAL(jcpoc-1, wp)
526         reminup = 1./ 400. * EXP(remindelta)
527         !
528         ! Discretization based on incomplete gamma functions
529         ! As incomplete gamma functions are not available in standard
530         ! fortran 95, they have been coded as functions in this module (gamain)
531         ! ---------------------------------------------------------------------
532         !
533         alphan(1) = gamain(reminup, rshape, ifault)
534         reminp(1) = gamain(reminup, rshape+1.0, ifault) * xremip / alphan(1)
535         DO jn = 2, jcpoc-1
536            reminup = 1./ 400. * EXP( REAL(jn, wp) * remindelta)
537            remindown = 1. / 400. * EXP( REAL(jn-1, wp) * remindelta)
538            alphan(jn) = gamain(reminup, rshape, ifault) - gamain(remindown, rshape, ifault)
539            reminp(jn) = gamain(reminup, rshape+1.0, ifault) - gamain(remindown, rshape+1.0, ifault)
540            reminp(jn) = reminp(jn) * xremip / alphan(jn)
541         END DO
542         remindown = 1. / 400. * EXP( REAL(jcpoc-1, wp) * remindelta)
543         alphan(jcpoc) = 1.0 - gamain(remindown, rshape, ifault)
544         reminp(jcpoc) = 1.0 - gamain(remindown, rshape+1.0, ifault)
545         reminp(jcpoc) = reminp(jcpoc) * xremip / alphan(jcpoc)
546
547      ELSE
548         alphan(jcpoc) = 1.
549         reminp(jcpoc) = xremip
550      ENDIF
551
552      DO jn = 1, jcpoc
553         alphap(:,:,:,jn) = alphan(jn)
554      END DO
555
556   END SUBROUTINE p4z_poc_init
557
558
559   REAL FUNCTION alngam( xvalue, ifault )
560      !*****************************************************************************80
561      !
562      !! ALNGAM computes the logarithm of the gamma function.
563      !
564      !  Modified:    13 January 2008
565      !
566      !  Author  :    Allan Macleod
567      !               FORTRAN90 version by John Burkardt
568      !
569      !  Reference:
570      !    Allan Macleod, Algorithm AS 245,
571      !    A Robust and Reliable Algorithm for the Logarithm of the Gamma Function,
572      !    Applied Statistics,
573      !    Volume 38, Number 2, 1989, pages 397-402.
574      !
575      !  Parameters:
576      !
577      !    Input, real ( kind = 8 ) XVALUE, the argument of the Gamma function.
578      !
579      !    Output, integer ( kind = 4 ) IFAULT, error flag.
580      !    0, no error occurred.
581      !    1, XVALUE is less than or equal to 0.
582      !    2, XVALUE is too big.
583      !
584      !    Output, real ( kind = 8 ) ALNGAM, the logarithm of the gamma function of X.
585      !*****************************************************************************80
586  implicit none
587
588  real(wp), parameter :: alr2pi = 0.918938533204673E+00
589  integer:: ifault
590  real(wp), dimension ( 9 ) :: r1 = (/ &
591    -2.66685511495E+00, &
592    -24.4387534237E+00, &
593    -21.9698958928E+00, &
594     11.1667541262E+00, &
595     3.13060547623E+00, &
596     0.607771387771E+00, &
597     11.9400905721E+00, &
598     31.4690115749E+00, &
599     15.2346874070E+00 /)
600  real(wp), dimension ( 9 ) :: r2 = (/ &
601    -78.3359299449E+00, &
602    -142.046296688E+00, &
603     137.519416416E+00, &
604     78.6994924154E+00, &
605     4.16438922228E+00, &
606     47.0668766060E+00, &
607     313.399215894E+00, &
608     263.505074721E+00, &
609     43.3400022514E+00 /)
610  real(wp), dimension ( 9 ) :: r3 = (/ &
611    -2.12159572323E+05, &
612     2.30661510616E+05, &
613     2.74647644705E+04, &
614    -4.02621119975E+04, &
615    -2.29660729780E+03, &
616    -1.16328495004E+05, &
617    -1.46025937511E+05, &
618    -2.42357409629E+04, &
619    -5.70691009324E+02 /)
620  real(wp), dimension ( 5 ) :: r4 = (/ &
621     0.279195317918525E+00, &
622     0.4917317610505968E+00, &
623     0.0692910599291889E+00, &
624     3.350343815022304E+00, &
625     6.012459259764103E+00 /)
626  real (wp) :: x
627  real (wp) :: x1
628  real (wp) :: x2
629  real (wp), parameter :: xlge = 5.10E+05
630  real (wp), parameter :: xlgst = 1.0E+30
631  real (wp) :: xvalue
632  real (wp) :: y
633
634  x = xvalue
635  alngam = 0.0E+00
636!
637!  Check the input.
638!
639  if ( xlgst <= x ) then
640    ifault = 2
641    return
642  end if
643  if ( x <= 0.0E+00 ) then
644    ifault = 1
645    return
646  end if
647
648  ifault = 0
649!
650!  Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined.
651!
652  if ( x < 1.5E+00 ) then
653
654    if ( x < 0.5E+00 ) then
655      alngam = - log ( x )
656      y = x + 1.0E+00
657!
658!  Test whether X < machine epsilon.
659!
660      if ( y == 1.0E+00 ) then
661        return
662      end if
663
664    else
665
666      alngam = 0.0E+00
667      y = x
668      x = ( x - 0.5E+00 ) - 0.5E+00
669
670    end if
671
672    alngam = alngam + x * (((( &
673        r1(5)   * y &
674      + r1(4) ) * y &
675      + r1(3) ) * y &
676      + r1(2) ) * y &
677      + r1(1) ) / (((( &
678                  y &
679      + r1(9) ) * y &
680      + r1(8) ) * y &
681      + r1(7) ) * y &
682      + r1(6) )
683
684    return
685
686  end if
687!
688!  Calculation for 1.5 <= X < 4.0.
689!
690  if ( x < 4.0E+00 ) then
691
692    y = ( x - 1.0E+00 ) - 1.0E+00
693
694    alngam = y * (((( &
695        r2(5)   * x &
696      + r2(4) ) * x &
697      + r2(3) ) * x &
698      + r2(2) ) * x &
699      + r2(1) ) / (((( &
700                  x &
701      + r2(9) ) * x &
702      + r2(8) ) * x &
703      + r2(7) ) * x &
704      + r2(6) )
705!
706!  Calculation for 4.0 <= X < 12.0.
707!
708  else if ( x < 12.0E+00 ) then
709
710    alngam = (((( &
711        r3(5)   * x &
712      + r3(4) ) * x &
713      + r3(3) ) * x &
714      + r3(2) ) * x &
715      + r3(1) ) / (((( &
716                  x &
717      + r3(9) ) * x &
718      + r3(8) ) * x &
719      + r3(7) ) * x &
720      + r3(6) )
721!
722!  Calculation for 12.0 <= X.
723!
724  else
725
726    y = log ( x )
727    alngam = x * ( y - 1.0E+00 ) - 0.5E+00 * y + alr2pi
728
729    if ( x <= xlge ) then
730
731      x1 = 1.0E+00 / x
732      x2 = x1 * x1
733
734      alngam = alngam + x1 * ( ( &
735             r4(3)   * &
736        x2 + r4(2) ) * &
737        x2 + r4(1) ) / ( ( &
738        x2 + r4(5) ) * &
739        x2 + r4(4) )
740
741    end if
742
743   end if
744
745   END FUNCTION alngam
746
747
748   REAL FUNCTION gamain( x, p, ifault )
749!*****************************************************************************80
750!
751!! GAMAIN computes the incomplete gamma ratio.
752!
753!  Discussion:
754!
755!    A series expansion is used if P > X or X <= 1.  Otherwise, a
756!    continued fraction approximation is used.
757!
758!  Modified:
759!
760!    17 January 2008
761!
762!  Author:
763!
764!    G Bhattacharjee
765!    FORTRAN90 version by John Burkardt
766!
767!  Reference:
768!
769!    G Bhattacharjee,
770!    Algorithm AS 32:
771!    The Incomplete Gamma Integral,
772!    Applied Statistics,
773!    Volume 19, Number 3, 1970, pages 285-287.
774!
775!  Parameters:
776!
777!    Input, real ( kind = 8 ) X, P, the parameters of the incomplete
778!    gamma ratio.  0 <= X, and 0 < P.
779!
780!    Output, integer ( kind = 4 ) IFAULT, error flag.
781!    0, no errors.
782!    1, P <= 0.
783!    2, X < 0.
784!    3, underflow.
785!    4, error return from the Log Gamma routine.
786!
787!    Output, real ( kind = 8 ) GAMAIN, the value of the incomplete
788!    gamma ratio.
789!
790  implicit none
791
792  real (wp) a
793  real (wp), parameter :: acu = 1.0E-08
794  real (wp) an
795  real (wp) arg
796  real (wp) b
797  real (wp) dif
798  real (wp) factor
799  real (wp) g
800  real (wp) gin
801  integer i
802  integer ifault
803  real (wp), parameter :: oflo = 1.0E+37
804  real (wp) p
805  real (wp) pn(6)
806  real (wp) rn
807  real (wp) term
808  real (wp), parameter :: uflo = 1.0E-37
809  real (wp) x
810!
811!  Check the input.
812!
813  if ( p <= 0.0E+00 ) then
814    ifault = 1
815    gamain = 0.0E+00
816    return
817  end if
818
819  if ( x < 0.0E+00 ) then
820    ifault = 2
821    gamain = 0.0E+00
822    return
823  end if
824
825  if ( x == 0.0E+00 ) then
826    ifault = 0
827    gamain = 0.0E+00
828    return
829  end if
830
831  g = alngam ( p, ifault )
832
833  if ( ifault /= 0 ) then
834    ifault = 4
835    gamain = 0.0E+00
836    return
837  end if
838
839  arg = p * log ( x ) - x - g
840
841  if ( arg < log ( uflo ) ) then
842    ifault = 3
843    gamain = 0.0E+00
844    return
845  end if
846
847  ifault = 0
848  factor = exp ( arg )
849!
850!  Calculation by series expansion.
851!
852  if ( x <= 1.0E+00 .or. x < p ) then
853
854    gin = 1.0E+00
855    term = 1.0E+00
856    rn = p
857
858    do
859
860      rn = rn + 1.0E+00
861      term = term * x / rn
862      gin = gin + term
863
864      if ( term <= acu ) then
865        exit
866      end if
867
868    end do
869
870    gamain = gin * factor / p
871    return
872
873  end if
874!
875!  Calculation by continued fraction.
876!
877  a = 1.0E+00 - p
878  b = a + x + 1.0E+00
879  term = 0.0E+00
880
881  pn(1) = 1.0E+00
882  pn(2) = x
883  pn(3) = x + 1.0E+00
884  pn(4) = x * b
885
886  gin = pn(3) / pn(4)
887
888  do
889
890    a = a + 1.0E+00
891    b = b + 2.0E+00
892    term = term + 1.0E+00
893    an = a * term
894    do i = 1, 2
895      pn(i+4) = b * pn(i+2) - an * pn(i)
896    end do
897
898    if ( pn(6) /= 0.0E+00 ) then
899
900      rn = pn(5) / pn(6)
901      dif = abs ( gin - rn )
902!
903!  Absolute error tolerance satisfied?
904!
905      if ( dif <= acu ) then
906!
907!  Relative error tolerance satisfied?
908!
909        if ( dif <= acu * rn ) then
910          gamain = 1.0E+00 - factor * gin
911          exit
912        end if
913
914      end if
915
916      gin = rn
917
918    end if
919
920    do i = 1, 4
921      pn(i) = pn(i+2)
922    end do
923    if ( oflo <= abs ( pn(5) ) ) then
924
925      do i = 1, 4
926        pn(i) = pn(i) / oflo
927      end do
928
929    end if
930
931  end do
932
933  END FUNCTION gamain
934
935   !!======================================================================
936END MODULE p4zpoc
Note: See TracBrowser for help on using the repository browser.