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 branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90 @ 9450

Last change on this file since 9450 was 9450, checked in by aumont, 6 years ago

debug PISCES code

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