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.
p4zsink.F90 in branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago

changes after review
File size: 48.1 KB
Line 
1MODULE p4zsink
2   !!======================================================================
3   !!                         ***  MODULE p4zsink  ***
4   !! TOP :  PISCES  vertical flux of particulate matter due to gravitational sinking
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) Change aggregation formula
9   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
10   !!----------------------------------------------------------------------
11#if defined key_pisces
12   !!----------------------------------------------------------------------
13   !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
14   !!   p4z_sink_init  :  Unitialisation of sinking speed parameters
15   !!   p4z_sink_alloc :  Allocate sinking speed variables
16   !!----------------------------------------------------------------------
17   USE oce_trc         !  shared variables between ocean and passive tracers
18   USE trc             !  passive tracers common variables
19   USE sms_pisces      !  PISCES Source Minus Sink variables
20   USE prtctl_trc      !  print control for debugging
21   USE iom             !  I/O manager
22   USE lib_mpp
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   p4z_sink         ! called in p4zbio.F90
28   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
29   PUBLIC   p4z_sink_alloc
30
31   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio3   !: POC sinking speed
32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio4   !: GOC sinking speed
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wscal    !: Calcite and BSi sinking speeds
34
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinking, sinking2  !: POC sinking fluxes
36   !                                                          !  (different meanings depending on the parameterization)
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkcal, sinksil   !: CaCO3 and BSi sinking fluxes
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
39#if ! defined key_kriest
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
41#endif
42
43   INTEGER  :: ik100
44
45#if  defined key_kriest
46   REAL(wp) ::  xkr_sfact    !: Sinking factor
47   REAL(wp) ::  xkr_stick    !: Stickiness
48   REAL(wp) ::  xkr_nnano    !: Nbr of cell in nano size class
49   REAL(wp) ::  xkr_ndiat    !: Nbr of cell in diatoms size class
50   REAL(wp) ::  xkr_nmicro   !: Nbr of cell in microzoo size class
51   REAL(wp) ::  xkr_nmeso    !: Nbr of cell in mesozoo  size class
52   REAL(wp) ::  xkr_naggr    !: Nbr of cell in aggregates  size class
53
54   REAL(wp) ::  xkr_frac 
55
56   REAL(wp), PUBLIC ::  xkr_dnano       !: Size of particles in nano pool
57   REAL(wp), PUBLIC ::  xkr_ddiat       !: Size of particles in diatoms pool
58   REAL(wp), PUBLIC ::  xkr_dmicro      !: Size of particles in microzoo pool
59   REAL(wp), PUBLIC ::  xkr_dmeso       !: Size of particles in mesozoo pool
60   REAL(wp), PUBLIC ::  xkr_daggr       !: Size of particles in aggregates pool
61   REAL(wp), PUBLIC ::  xkr_wsbio_min   !: min vertical particle speed
62   REAL(wp), PUBLIC ::  xkr_wsbio_max   !: max vertical particle speed
63
64   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   xnumm   !:  maximum number of particles in aggregates
65#endif
66
67   !!----------------------------------------------------------------------
68   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
69   !! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74#if ! defined key_kriest
75   !!----------------------------------------------------------------------
76   !!   'standard sinking parameterisation'                  ???
77   !!----------------------------------------------------------------------
78
79   SUBROUTINE p4z_sink ( kt, knt )
80      !!---------------------------------------------------------------------
81      !!                     ***  ROUTINE p4z_sink  ***
82      !!
83      !! ** Purpose :   Compute vertical flux of particulate matter due to
84      !!                gravitational sinking
85      !!
86      !! ** Method  : - ???
87      !!---------------------------------------------------------------------
88      INTEGER, INTENT(in) :: kt, knt
89      INTEGER  ::   ji, jj, jk, jit
90      INTEGER  ::   iiter1, iiter2
91      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
92      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
93      REAL(wp) ::   zfact, zwsmax, zmax, zstep
94      CHARACTER (len=25) :: charout
95      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
96      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
97      !!---------------------------------------------------------------------
98      !
99      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
100      !
101      !    Sinking speeds of detritus is increased with depth as shown
102      !    by data and from the coagulation theory
103      !    -----------------------------------------------------------
104!$OMP PARALLEL
105!$OMP DO schedule(static) private(jk, jj, ji, zmax, zfact)
106      DO jk = 1, jpkm1
107         DO jj = 1, jpj
108            DO ji = 1,jpi
109               zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
110               zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / 5000._wp
111               wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
112            END DO
113         END DO
114      END DO
115
116      ! limit the values of the sinking speeds to avoid numerical instabilities 
117!$OMP DO schedule(static) private(jk, jj, ji)
118      DO jk = 1, jpk
119         DO jj = 1, jpj
120            DO ji = 1, jpi
121               wsbio3(ji,jj,jk) = wsbio
122               wscal (ji,jj,jk) = wsbio4(ji,jj,jk)
123            END DO
124         END DO
125      END DO
126!$OMP END PARALLEL
127      !
128      ! OA This is (I hope) a temporary solution for the problem that may
129      ! OA arise in specific situation where the CFL criterion is broken
130      ! OA for vertical sedimentation of particles. To avoid this, a time
131      ! OA splitting algorithm has been coded. A specific maximum
132      ! OA iteration number is provided and may be specified in the namelist
133      ! OA This is to avoid very large iteration number when explicit free
134      ! OA surface is used (for instance). When niter?max is set to 1,
135      ! OA this computation is skipped. The crude old threshold method is
136      ! OA then applied. This also happens when niter exceeds nitermax.
137      IF( MAX( niter1max, niter2max ) == 1 ) THEN
138        iiter1 = 1
139        iiter2 = 1
140      ELSE
141        iiter1 = 1
142        iiter2 = 1
143!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zwsmax, iiter1, iiter2)
144        DO jk = 1, jpkm1
145          DO jj = 1, jpj
146             DO ji = 1, jpi
147                IF( tmask(ji,jj,jk) == 1) THEN
148                   zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep
149                   iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
150                   iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
151                ENDIF
152             END DO
153          END DO
154        END DO
155        IF( lk_mpp ) THEN
156           CALL mpp_max( iiter1 )
157           CALL mpp_max( iiter2 )
158        ENDIF
159        iiter1 = MIN( iiter1, niter1max )
160        iiter2 = MIN( iiter2, niter2max )
161      ENDIF
162
163!$OMP PARALLEL
164!$OMP DO schedule(static) private(jk, jj, ji, zwsmax)
165      DO jk = 1,jpkm1
166         DO jj = 1, jpj
167            DO ji = 1, jpi
168               IF( tmask(ji,jj,jk) == 1 ) THEN
169                 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
170                 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
171                 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
172               ENDIF
173            END DO
174         END DO
175      END DO
176!$OMP END DO NOWAIT
177
178      !  Initializa to zero all the sinking arrays
179      !   -----------------------------------------
180!$OMP DO schedule(static) private(jk, jj, ji)
181      DO jk = 1, jpk
182         DO jj = 1, jpj
183            DO ji = 1, jpi
184               sinking (ji,jj,jk) = 0.e0
185               sinking2(ji,jj,jk) = 0.e0
186               sinkcal (ji,jj,jk) = 0.e0
187               sinkfer (ji,jj,jk) = 0.e0
188               sinksil (ji,jj,jk) = 0.e0
189               sinkfer2(ji,jj,jk) = 0.e0
190            END DO
191         END DO
192      END DO
193!$OMP END PARALLEL
194
195      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
196      !   -----------------------------------------------------
197      DO jit = 1, iiter1
198        CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
199        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
200      END DO
201
202      DO jit = 1, iiter2
203        CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
204        CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
205        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
206        CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
207      END DO
208
209      !  Exchange between organic matter compartments due to coagulation/disaggregation
210      !  ---------------------------------------------------
211!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zstep, zfact, zagg, zagg1, zagg2) &
212!$OMP& private(zagg3, zagg4, zaggfe, zaggdoc, zaggdoc2, zaggdoc3)
213      DO jk = 1, jpkm1
214         DO jj = 1, jpj
215            DO ji = 1, jpi
216               !
217               zstep = xstep 
218# if defined key_degrad
219               zstep = zstep * facvol(ji,jj,jk)
220# endif
221               zfact = zstep * xdiss(ji,jj,jk)
222               !  Part I : Coagulation dependent on turbulence
223               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
224               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
225
226               ! Part II : Differential settling
227
228               !  Aggregation of small into large particles
229               zagg3 =  47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
230               zagg4 =  3.3  * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
231
232               zagg   = zagg1 + zagg2 + zagg3 + zagg4
233               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
234
235               ! Aggregation of DOC to POC :
236               ! 1st term is shear aggregation of DOC-DOC
237               ! 2nd term is shear aggregation of DOC-POC
238               ! 3rd term is differential settling of DOC-POC
239               zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
240               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
241               ! transfer of DOC to GOC :
242               ! 1st term is shear aggregation
243               ! 2nd term is differential settling
244               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
245               ! tranfer of DOC to POC due to brownian motion
246               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
247
248               !  Update the trends
249               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
250               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
251               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
252               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
253               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
254               !
255            END DO
256         END DO
257      END DO
258
259
260     ! Total carbon export per year
261     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  ) &
262        &        t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
263     !
264     IF( lk_iomput ) THEN
265       IF( knt == nrdttrc ) THEN
266          CALL wrk_alloc( jpi, jpj,      zw2d )
267          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
268          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
269          !
270          IF( iom_use( "EPC100" ) )  THEN
271!$OMP DO schedule(static) private(jj, ji)
272             DO jj = 1, jpj
273                DO ji = 1, jpi
274                   zw2d(ji,jj) = ( sinking(ji,jj,ik100) + sinking2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
275                END DO
276             END DO
277             CALL iom_put( "EPC100"  , zw2d )
278          ENDIF
279          IF( iom_use( "EPFE100" ) )  THEN
280!$OMP DO schedule(static) private(jj, ji)
281             DO jj = 1, jpj
282                DO ji = 1, jpi
283                   zw2d(ji,jj) = ( sinkfer(ji,jj,ik100) + sinkfer2(ji,jj,ik100) ) * zfact * tmask(ji,jj,1) ! Export of iron at 100m
284                END DO
285             END DO
286             CALL iom_put( "EPFE100"  , zw2d )
287          ENDIF
288          IF( iom_use( "EPCAL100" ) )  THEN
289!$OMP DO schedule(static) private(jj, ji)
290             DO jj = 1, jpj
291                DO ji = 1, jpi
292                   zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of calcite at 100m
293                END DO
294             END DO
295             CALL iom_put( "EPCAL100"  , zw2d )
296          ENDIF
297          IF( iom_use( "EPSI100" ) )  THEN
298!$OMP DO schedule(static) private(jj, ji)
299             DO jj = 1, jpj
300                DO ji = 1, jpi
301                   zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
302                END DO
303             END DO
304             CALL iom_put( "EPSI100"  , zw2d )
305          ENDIF
306          IF( iom_use( "EXPC" ) )  THEN
307!$OMP DO schedule(static) private(jk, jj, ji)
308             DO jk = 1, jpk
309                DO jj = 1, jpj
310                   DO ji = 1, jpi
311                      zw3d(ji,jj,jk) = ( sinking(ji,jj,jk) + sinking2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
312                   END DO
313                END DO
314             END DO
315             CALL iom_put( "EXPC"  , zw3d )
316          ENDIF
317          IF( iom_use( "EXPFE" ) )  THEN
318!$OMP DO schedule(static) private(jk, jj, ji)
319             DO jk = 1, jpk
320                DO jj = 1, jpj
321                   DO ji = 1, jpi
322                      zw3d(ji,jj,jk) = ( sinkfer(ji,jj,jk) + sinkfer2(ji,jj,jk) ) * zfact * tmask(ji,jj,jk) ! Export of iron
323                   END DO
324                END DO
325             END DO
326             CALL iom_put( "EXPFE"  , zw3d )
327          ENDIF
328          IF( iom_use( "EXPCAL" ) )  THEN
329!$OMP DO schedule(static) private(jk, jj, ji)
330             DO jk = 1, jpk
331                DO jj = 1, jpj
332                   DO ji = 1, jpi
333                      zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite
334                   END DO
335                END DO
336             END DO
337             CALL iom_put( "EXPCAL"  , zw3d )
338          ENDIF
339          IF( iom_use( "EXPSI" ) )  THEN
340!$OMP DO schedule(static) private(jk, jj, ji)
341             DO jk = 1, jpk
342                DO jj = 1, jpj
343                   DO ji = 1, jpi
344                      zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
345                   END DO
346                END DO
347             END DO
348             CALL iom_put( "EXPSI"  , zw3d )
349          ENDIF
350          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
351          !
352          CALL wrk_dealloc( jpi, jpj,      zw2d )
353          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
354        ENDIF
355      ELSE
356         IF( ln_diatrc ) THEN
357            zfact = 1.e3 * rfact2r
358!$OMP DO schedule(static) private(jj, ji)
359            DO jj = 1, jpj
360               DO ji = 1, jpi
361                  trc2d(ji,jj,jp_pcs0_2d + 4) = sinking (ji,jj,ik100) * zfact * tmask(ji,jj,1)
362                  trc2d(ji,jj,jp_pcs0_2d + 5) = sinking2(ji,jj,ik100) * zfact * tmask(ji,jj,1)
363                  trc2d(ji,jj,jp_pcs0_2d + 6) = sinkfer (ji,jj,ik100) * zfact * tmask(ji,jj,1)
364                  trc2d(ji,jj,jp_pcs0_2d + 7) = sinkfer2(ji,jj,ik100) * zfact * tmask(ji,jj,1)
365                  trc2d(ji,jj,jp_pcs0_2d + 8) = sinksil (ji,jj,ik100) * zfact * tmask(ji,jj,1)
366                  trc2d(ji,jj,jp_pcs0_2d + 9) = sinkcal (ji,jj,ik100) * zfact * tmask(ji,jj,1)
367               END DO
368            END DO
369         ENDIF
370      ENDIF
371      !
372      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
373         WRITE(charout, FMT="('sink')")
374         CALL prt_ctl_trc_info(charout)
375         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
376      ENDIF
377      !
378      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
379      !
380   END SUBROUTINE p4z_sink
381
382   SUBROUTINE p4z_sink_init
383      !!----------------------------------------------------------------------
384      !!                  ***  ROUTINE p4z_sink_init  ***
385      !!----------------------------------------------------------------------
386      INTEGER :: jk
387
388      ik100 = 10        !  last level where depth less than 100 m
389      DO jk = jpkm1, 1, -1
390         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
391      END DO
392      IF (lwp) WRITE(numout,*)
393      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
394      IF (lwp) WRITE(numout,*)
395      !
396      t_oce_co2_exp = 0._wp
397      !
398   END SUBROUTINE p4z_sink_init
399
400#else
401   !!----------------------------------------------------------------------
402   !!   'Kriest sinking parameterisation'        key_kriest          ???
403   !!----------------------------------------------------------------------
404
405   SUBROUTINE p4z_sink ( kt, knt )
406      !!---------------------------------------------------------------------
407      !!                ***  ROUTINE p4z_sink  ***
408      !!
409      !! ** Purpose :   Compute vertical flux of particulate matter due to
410      !!              gravitational sinking - Kriest parameterization
411      !!
412      !! ** Method  : - ???
413      !!---------------------------------------------------------------------
414      !
415      INTEGER, INTENT(in) :: kt, knt
416      !
417      INTEGER  :: ji, jj, jk, jit, niter1, niter2
418      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh
419      REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc
420      REAL(wp) :: znum , zeps, zfm, zgm, zsm
421      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
422      REAL(wp) :: zval1, zval2, zval3, zval4
423      REAL(wp) :: zfact
424      INTEGER  :: ik1
425      CHARACTER (len=25) :: charout
426      REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d 
427      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
428      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
429      !!---------------------------------------------------------------------
430      !
431      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
432      !
433      CALL wrk_alloc( jpi, jpj, jpk, znum3d )
434      !
435      !     Initialisation of variables used to compute Sinking Speed
436      !     ---------------------------------------------------------
437
438      zval1 = 1. + xkr_zeta
439      zval2 = 1. + xkr_zeta + xkr_eta
440      zval3 = 1. + xkr_eta
441!$OMP PARALLEL
442!$OMP DO schedule(static) private(jk, jj, ji)
443      DO jk = 1, jpk
444         DO jj = 1, jpj
445            DO ji = 1, jpi
446               znum3d(ji,jj,jk) = 0.e0
447            END DO
448         END DO
449      END DO
450      !     Computation of the vertical sinking speed : Kriest et Evans, 2000
451      !     -----------------------------------------------------------------
452
453!$OMP DO schedule(static) private(jk, jj, ji, znum, zeps, zfm, zgm, zdiv, zdiv1)
454      DO jk = 1, jpkm1
455         DO jj = 1, jpj
456            DO ji = 1, jpi
457               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
458                  znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
459                  ! -------------- To avoid sinking speed over 50 m/day -------
460                  znum  = MIN( xnumm(jk), znum )
461                  znum  = MAX( 1.1      , znum )
462                  znum3d(ji,jj,jk) = znum
463                  !------------------------------------------------------------
464                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
465                  zfm   = xkr_frac**( 1. - zeps )
466                  zgm   = xkr_frac**( zval1 - zeps )
467                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
468                  zdiv1 = zeps - zval3
469                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
470                     &             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
471                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
472                     &             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
473                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
474               ENDIF
475            END DO
476         END DO
477      END DO
478!$OMP END DO NOWAIT
479!$OMP DO schedule(static) private(jk, jj, ji)
480      DO jk = 1, jpk
481         DO jj = 1, jpj
482            DO ji = 1, jpi
483               wscal(ji,jj,jk) = MAX( wsbio3(ji,jj,jk), 30._wp )
484            END DO
485         END DO
486      END DO
487!$OMP END DO NOWAIT
488
489      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
490      !   -----------------------------------------
491!$OMP DO schedule(static) private(jk, jj, ji)
492      DO jk = 1, jpk
493         DO jj = 1, jpj
494            DO ji = 1, jpi
495               sinking (ji,jj,jk) = 0.e0
496               sinking2(ji,jj,jk) = 0.e0
497               sinkcal (ji,jj,jk) = 0.e0
498               sinkfer (ji,jj,jk) = 0.e0
499               sinksil (ji,jj,jk) = 0.e0
500            END DO
501         END DO
502      END DO
503!$OMP END PARALLEL
504     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
505     !   -----------------------------------------------------
506
507      niter1 = niter1max
508      niter2 = niter2max
509
510      DO jit = 1, niter1
511        CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )
512        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )
513        CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
514        CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
515      END DO
516
517      DO jit = 1, niter2
518        CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )
519      END DO
520
521     !  Exchange between organic matter compartments due to coagulation/disaggregation
522     !  ---------------------------------------------------
523
524      zval1 = 1. + xkr_zeta
525      zval2 = 1. + xkr_eta
526      zval3 = 3. + xkr_eta
527      zval4 = 4. + xkr_eta
528
529!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, znum, zeps, zfm, zsm, zdiv, zdiv1, zdiv2, zdiv3, zdiv4) &
530!$OMP& private(zdiv5, zagg, zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggdoc, zaggdoc1, zaggsh, zaggsi, znumdoc)
531      DO jk = 1,jpkm1
532         DO jj = 1,jpj
533            DO ji = 1,jpi
534               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
535
536                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp
537                  !-------------- To avoid sinking speed over 50 m/day -------
538                  znum  = min(xnumm(jk),znum)
539                  znum  = MAX( 1.1,znum)
540                  !------------------------------------------------------------
541                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
542                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
543                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
544                  zdiv2 = zeps - 2.
545                  zdiv3 = zeps - 3.
546                  zdiv4 = zeps - zval2
547                  zdiv5 = 2.* zeps - zval4
548                  zfm   = xkr_frac**( 1.- zeps )
549                  zsm   = xkr_frac**xkr_eta
550
551                  !    Part I : Coagulation dependant on turbulence
552                  !    ----------------------------------------------
553
554                  zagg1 =  0.163 * trb(ji,jj,jk,jpnum)**2               &
555                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
556                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
557                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
558                     &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
559                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm*                       &
560                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
561                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
562                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
563                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
564                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
565                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))   
566
567                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3 
568                 
569                 !    Aggregation of small into large particles
570                 !    Part II : Differential settling
571                 !    ----------------------------------------------
572
573                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       &
574                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
575                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
576                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
577                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
578                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
579
580                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         &
581                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
582                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
583                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
584                     &                 /zdiv) 
585
586                  !
587                  !     Fractionnation by swimming organisms
588                  !     ------------------------------------
589
590                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  &
591                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
592                    &      * 10000.*xstep
593
594                  !     Aggregation of DOC to small particles
595                  !     --------------------------------------
596
597                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   &
598                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)
599                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  &
600                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)
601
602# if defined key_degrad
603                   zagg1   = zagg1   * facvol(ji,jj,jk)                 
604                   zagg2   = zagg2   * facvol(ji,jj,jk)                 
605                   zagg3   = zagg3   * facvol(ji,jj,jk)                 
606                   zagg4   = zagg4   * facvol(ji,jj,jk)                 
607                   zagg5   = zagg5   * facvol(ji,jj,jk)                 
608                   zaggdoc = zaggdoc * facvol(ji,jj,jk)                 
609                   zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)
610# endif
611                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
612                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
613                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
614                  !
615                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
616                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
617                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
618                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1
619
620               ENDIF
621            END DO
622         END DO
623      END DO
624
625     ! Total primary production per year
626     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )
627     !
628     IF( lk_iomput ) THEN
629        IF( knt == nrdttrc ) THEN
630          CALL wrk_alloc( jpi, jpj,      zw2d )
631          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
632          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
633          !
634          IF( iom_use( "EPC100" ) )  THEN
635!$OMP DO schedule(static) private(jj, ji)
636             DO jj = 1, jpj
637                DO ji = 1, jpi
638                   zw2d(ji,jj) = sinking(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of carbon at 100m
639                END DO
640             END DO
641             CALL iom_put( "EPC100"  , zw2d )
642          ENDIF
643          IF( iom_use( "EPN100" ) )  THEN
644!$OMP DO schedule(static) private(jj, ji)
645             DO jj = 1, jpj
646                DO ji = 1, jpi
647                   zw2d(ji,jj) = sinking2(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of number of aggregates ?
648                END DO
649             END DO
650             CALL iom_put( "EPN100"  , zw2d )
651          ENDIF
652          IF( iom_use( "EPCAL100" ) )  THEN
653!$OMP DO schedule(static) private(jj, ji)
654             DO jj = 1, jpj
655                DO ji = 1, jpi
656                   zw2d(ji,jj) = sinkcal(ji,jj,ik100) * zfact * tmask(ji,jj,1) !Export of calcite at 100m
657                END DO
658             END DO
659             CALL iom_put( "EPCAL100"  , zw2d )
660          ENDIF
661          IF( iom_use( "EPSI100" ) )  THEN
662!$OMP DO schedule(static) private(jj, ji)
663             DO jj = 1, jpj
664                DO ji = 1, jpi
665                   zw2d(ji,jj) =  sinksil(ji,jj,ik100) * zfact * tmask(ji,jj,1) ! Export of bigenic silica at 100m
666                END DO
667             END DO
668             CALL iom_put( "EPSI100"  , zw2d )
669          ENDIF
670          IF( iom_use( "EXPC" ) )  THEN
671!$OMP DO schedule(static) private(jk, jj, ji)
672             DO jk = 1, jpk
673                DO jj = 1, jpj
674                   DO ji = 1, jpi
675                      zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
676                   END DO
677                END DO
678             END DO
679             CALL iom_put( "EXPC"  , zw3d )
680          ENDIF
681          IF( iom_use( "EXPN" ) )  THEN
682!$OMP DO schedule(static) private(jk, jj, ji)
683             DO jk = 1, jpk
684                DO jj = 1, jpj
685                   DO ji = 1, jpi
686                      zw3d(ji,jj,jk) = sinking(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of carbon in the water column
687                   END DO
688                END DO
689             END DO
690             CALL iom_put( "EXPN"  , zw3d )
691          ENDIF
692          IF( iom_use( "EXPCAL" ) )  THEN
693!$OMP DO schedule(static) private(jk, jj, ji)
694             DO jk = 1, jpk
695                DO jj = 1, jpj
696                   DO ji = 1, jpi
697                      zw3d(ji,jj,jk) = sinkcal(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of calcite
698                   END DO
699                END DO
700             END DO
701             CALL iom_put( "EXPCAL"  , zw3d )
702          ENDIF
703          IF( iom_use( "EXPSI" ) )  THEN
704!$OMP DO schedule(static) private(jk, jj, ji)
705             DO jk = 1, jpk
706                DO jj = 1, jpj
707                   DO ji = 1, jpi
708                      zw3d(ji,jj,jk) = sinksil(ji,jj,jk) * zfact * tmask(ji,jj,jk) ! Export of bigenic silica
709                   END DO
710                END DO
711             END DO
712             CALL iom_put( "EXPSI"  , zw3d )
713          ENDIF
714          IF( iom_use( "XNUM" ) )  THEN
715!$OMP DO schedule(static) private(jk, jj, ji)
716             DO jk = 1, jpk
717                DO jj = 1, jpj
718                   DO ji = 1, jpi
719                      zw3d(ji,jj,jk) = znum3d(ji,jj,jk) * tmask(ji,jj,jk) !  Number of particles on aggregats
720                   END DO
721                END DO
722             END DO
723             CALL iom_put( "XNUM"  , zw3d )
724          ENDIF
725          IF( iom_use( "WSC" ) )  THEN
726!$OMP DO schedule(static) private(jk, jj, ji)
727             DO jk = 1, jpk
728                DO jj = 1, jpj
729                   DO ji = 1, jpi
730                      zw3d(ji,jj,jk) = wsbio3(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of carbon particles
731                   END DO
732                END DO
733             END DO
734             CALL iom_put( "WSC"  , zw3d )
735          ENDIF
736          IF( iom_use( "WSN" ) )  THEN
737!$OMP DO schedule(static) private(jk, jj, ji)
738             DO jk = 1, jpk
739                DO jj = 1, jpj
740                   DO ji = 1, jpi
741                      zw3d(ji,jj,jk) = wsbio4(ji,jj,jk) * tmask(ji,jj,jk) ! Sinking speed of particles number
742                   END DO
743                END DO
744             END DO
745             CALL iom_put( "WSN"  , zw3d )
746          ENDIF
747          !
748          CALL wrk_dealloc( jpi, jpj,      zw2d )
749          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
750      ELSE
751         IF( ln_diatrc ) THEN
752            zfact = 1.e3 * rfact2r
753!$OMP PARALLEL
754!$OMP DO schedule(static) private(jj, ji)
755            DO jj = 1, jpj
756               DO ji = 1, jpi
757                  trc2d(ji,jj ,jp_pcs0_2d + 4)  = sinking (ji,jj,ik100)  * zfact * tmask(ji,jj,1)
758                  trc2d(ji,jj ,jp_pcs0_2d + 5)  = sinking2(ji,jj,ik100)  * zfact * tmask(ji,jj,1)
759                  trc2d(ji,jj ,jp_pcs0_2d + 6)  = sinkfer (ji,jj,ik100)  * zfact * tmask(ji,jj,1)
760                  trc2d(ji,jj ,jp_pcs0_2d + 7)  = sinksil (ji,jj,ik100)  * zfact * tmask(ji,jj,1)
761                  trc2d(ji,jj ,jp_pcs0_2d + 8)  = sinkcal (ji,jj,ik100)  * zfact * tmask(ji,jj,1)
762               END DO
763            END DO
764!$OMP DO schedule(static) private(jk, jj, ji)
765            DO jk = 1, jpk
766               DO jj = 1, jpj
767                  DO ji = 1, jpi
768                     trc3d(ji,jj,jk,jp_pcs0_3d + 11) = sinking (ji,jj,jk)      * zfact * tmask(ji,jj,jk)
769                     trc3d(ji,jj,jk,jp_pcs0_3d + 12) = sinking2(ji,jj,jk)      * zfact * tmask(ji,jj,jk)
770                     trc3d(ji,jj,jk,jp_pcs0_3d + 13) = sinksil (ji,jj,jk)      * zfact * tmask(ji,jj,jk)
771                     trc3d(ji,jj,jk,jp_pcs0_3d + 14) = sinkcal (ji,jj,jk)      * zfact * tmask(ji,jj,jk)
772                     trc3d(ji,jj,jk,jp_pcs0_3d + 15) = znum3d  (ji,jj,jk)              * tmask(ji,jj,jk)
773                     trc3d(ji,jj,jk,jp_pcs0_3d + 16) = wsbio3  (ji,jj,jk)              * tmask(ji,jj,jk)
774                     trc3d(ji,jj,jk,jp_pcs0_3d + 17) = wsbio4  (ji,jj,jk)              * tmask(ji,jj,jk)
775                  END DO
776               END DO
777            END DO
778!$OMP END PARALLEL
779         ENDIF
780      ENDIF
781
782      !
783      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
784         WRITE(charout, FMT="('sink')")
785         CALL prt_ctl_trc_info(charout)
786         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
787      ENDIF
788      !
789      CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
790      !
791      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
792      !
793   END SUBROUTINE p4z_sink
794
795
796   SUBROUTINE p4z_sink_init
797      !!----------------------------------------------------------------------
798      !!                  ***  ROUTINE p4z_sink_init  ***
799      !!
800      !! ** Purpose :   Initialization of sinking parameters
801      !!                Kriest parameterization only
802      !!
803      !! ** Method  :   Read the nampiskrs namelist and check the parameters
804      !!      called at the first timestep
805      !!
806      !! ** input   :   Namelist nampiskrs
807      !!----------------------------------------------------------------------
808      INTEGER  ::   jk, jn, kiter
809      INTEGER  ::   ios                 ! Local integer output status for namelist read
810      REAL(wp) ::   znum, zdiv
811      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
812      REAL(wp) ::   zmin, zmax, zl, zr, xacc
813      !
814      NAMELIST/nampiskrs/ xkr_sfact, xkr_stick ,  &
815         &                xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr
816      !!----------------------------------------------------------------------
817      !
818      IF( nn_timing == 1 )  CALL timing_start('p4z_sink_init')
819      !
820
821      REWIND( numnatp_ref )              ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest
822      READ  ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)
823901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )
824
825      REWIND( numnatp_cfg )              ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest
826      READ  ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )
827902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )
828      IF(lwm) WRITE ( numonp, nampiskrs )
829
830      IF(lwp) THEN
831         WRITE(numout,*)
832         WRITE(numout,*) ' Namelist : nampiskrs'
833         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
834         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
835         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
836         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
837         WRITE(numout,*) '    Nbr of cell in microzoo size class       xkr_nmicro   = ', xkr_nmicro
838         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
839         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
840      ENDIF
841
842
843      ! max and min vertical particle speed
844      xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
845      xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
846      IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
847
848      !
849      !    effect of the sizes of the different living pools on particle numbers
850      !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
851      !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
852      !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
853      !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
854      !    doc aggregates = 1um
855      ! ----------------------------------------------------------
856
857      xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
858      xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
859      xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )
860      xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
861      xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
862
863      !!---------------------------------------------------------------------
864      !!    'key_kriest'                                                  ???
865      !!---------------------------------------------------------------------
866      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
867      !  Search of the maximum number of particles in aggregates for each k-level.
868      !  Bissection Method
869      !--------------------------------------------------------------------
870      IF (lwp) THEN
871        WRITE(numout,*)
872        WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
873      ENDIF
874
875      xacc     =  0.001_wp
876      kiter    = 50
877      zmin     =  1.10_wp
878      zmax     = xkr_mass_max / xkr_mass_min
879      xkr_frac = zmax
880
881      DO jk = 1,jpk
882         zl = zmin
883         zr = zmax
884         wmax = 0.5 * e3t_n(1,1,jk) * rday * float(niter1max) / rfact2
885         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
886         znum = zl - 1.
887         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
888            & - ( xkr_wsbio_max * xkr_eta * znum * &
889            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
890            & - wmax
891
892         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
893         znum = zr - 1.
894         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
895            & - ( xkr_wsbio_max * xkr_eta * znum * &
896            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
897            & - wmax
898iflag:   DO jn = 1, kiter
899            IF    ( zwl == 0._wp ) THEN   ;   znummax = zl
900            ELSEIF( zwr == 0._wp ) THEN   ;   znummax = zr
901            ELSE
902               znummax = ( zr + zl ) / 2.
903               zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
904               znum = znummax - 1.
905               zws =  xkr_wsbio_min * xkr_zeta / zdiv &
906                  & - ( xkr_wsbio_max * xkr_eta * znum * &
907                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
908                  & - wmax
909               IF( zws * zwl < 0. ) THEN   ;   zr = znummax
910               ELSE                        ;   zl = znummax
911               ENDIF
912               zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
913               znum = zl - 1.
914               zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
915                  & - ( xkr_wsbio_max * xkr_eta * znum * &
916                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
917                  & - wmax
918
919               zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
920               znum = zr - 1.
921               zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
922                  & - ( xkr_wsbio_max * xkr_eta * znum * &
923                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
924                  & - wmax
925               !
926               IF ( ABS ( zws )  <= xacc ) EXIT iflag
927               !
928            ENDIF
929            !
930         END DO iflag
931
932         xnumm(jk) = znummax
933         IF (lwp) WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
934         !
935      END DO
936      !
937      ik100 = 10        !  last level where depth less than 100 m
938      DO jk = jpkm1, 1, -1
939         IF( gdept_1d(jk) > 100. )  iksed = jk - 1
940      END DO
941      IF (lwp) WRITE(numout,*)
942      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
943      IF (lwp) WRITE(numout,*)
944      !
945      t_oce_co2_exp = 0._wp
946      !
947      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink_init')
948      !
949  END SUBROUTINE p4z_sink_init
950
951#endif
952
953   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
954      !!---------------------------------------------------------------------
955      !!                     ***  ROUTINE p4z_sink2  ***
956      !!
957      !! ** Purpose :   Compute the sedimentation terms for the various sinking
958      !!     particles. The scheme used to compute the trends is based
959      !!     on MUSCL.
960      !!
961      !! ** Method  : - this ROUTINE compute not exactly the advection but the
962      !!      transport term, i.e.  div(u*tra).
963      !!---------------------------------------------------------------------
964      !
965      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
966      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting
967      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
968      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
969      !!
970      INTEGER  ::   ji, jj, jk, jn
971      REAL(wp) ::   zigma,zew,zign, zflx, zstep
972      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb 
973      !!---------------------------------------------------------------------
974      !
975      IF( nn_timing == 1 )  CALL timing_start('p4z_sink2')
976      !
977      ! Allocate temporary workspace
978      CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
979
980      zstep = rfact2 / FLOAT( kiter ) / 2.
981!$OMP PARALLEL
982!$OMP DO schedule(static) private(jk, jj, ji)
983      DO jk = 1, jpk
984         DO jj = 1, jpj
985            DO ji = 1, jpi
986               ztraz(ji,jj,jk) = 0.e0
987               zakz (ji,jj,jk) = 0.e0
988               ztrb (ji,jj,jk) = trb(ji,jj,jk,jp_tra)
989            END DO
990         END DO
991      END DO
992!$OMP END DO NOWAIT
993!$OMP DO schedule(static) private(jk)
994      DO jk = 1, jpkm1
995         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
996      END DO
997
998!$OMP DO schedule(static) private(jj, ji)
999      DO jj = 1, jpj
1000         DO ji = 1, jpi
1001            zwsink2(ji,jj,1) = 0.e0
1002         END DO
1003      END DO
1004!$OMP END DO NOWAIT
1005!$OMP END PARALLEL
1006
1007      IF( lk_degrad ) THEN
1008!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
1009         DO jk = 1, jpk
1010            DO jj = 1, jpj
1011               DO ji = 1, jpi
1012                  zwsink2(ji,jj,jk) = zwsink2(ji,jj,jk) * facvol(ji,jj,jk)
1013               END DO
1014            END DO
1015         END DO
1016      ENDIF
1017
1018
1019      ! Vertical advective flux
1020      DO jn = 1, 2
1021         !  first guess of the slopes interior values
1022!$OMP PARALLEL
1023!$OMP DO schedule(static) private(jk)
1024         DO jk = 2, jpkm1
1025            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
1026         END DO
1027!$OMP END DO NOWAIT
1028
1029!$OMP DO schedule(static) private(jj, ji)
1030      DO jj = 1, jpj
1031         DO ji = 1, jpi
1032            ztraz(ji,jj,1  ) = 0.0
1033            ztraz(ji,jj,jpk) = 0.0
1034         END DO
1035      END DO
1036
1037         ! slopes
1038!$OMP DO schedule(static) private(jk, jj, ji, zign)
1039         DO jk = 2, jpkm1
1040            DO jj = 1,jpj
1041               DO ji = 1, jpi
1042                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
1043                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
1044               END DO
1045            END DO
1046         END DO
1047         
1048         ! Slopes limitation
1049!$OMP DO schedule(static) private(jk, jj, ji)
1050         DO jk = 2, jpkm1
1051            DO jj = 1, jpj
1052               DO ji = 1, jpi
1053                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
1054                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
1055               END DO
1056            END DO
1057         END DO
1058         
1059         ! vertical advective flux
1060!$OMP DO schedule(static) private(jk, jj, ji, zigma, zew)
1061         DO jk = 1, jpkm1
1062            DO jj = 1, jpj     
1063               DO ji = 1, jpi   
1064                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
1065                  zew   = zwsink2(ji,jj,jk+1)
1066                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
1067               END DO
1068            END DO
1069         END DO
1070         !
1071!$OMP DO schedule(static) private(jj, ji)
1072         DO jj = 1, jpj
1073            DO ji = 1, jpi
1074                ! Boundary conditions
1075               psinkflx(ji,jj,1  ) = 0.e0
1076               psinkflx(ji,jj,jpk) = 0.e0
1077            END DO
1078         END DO
1079         
1080!$OMP DO schedule(static) private(jk, jj, ji, zflx)
1081         DO jk=1,jpkm1
1082            DO jj = 1,jpj
1083               DO ji = 1, jpi
1084                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
1085                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
1086               END DO
1087            END DO
1088         END DO
1089!$OMP END DO NOWAIT
1090!$OMP END PARALLEL
1091      ENDDO
1092!$OMP PARALLEL
1093!$OMP DO schedule(static) private(jk, jj, ji, zflx)
1094      DO jk = 1,jpkm1
1095         DO jj = 1,jpj
1096            DO ji = 1, jpi
1097               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
1098               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
1099            END DO
1100         END DO
1101      END DO
1102
1103!$OMP DO schedule(static) private(jk, jj, ji)
1104      DO jk = 1, jpk
1105         DO jj = 1, jpj
1106            DO ji = 1, jpi
1107               trb(ji,jj,jk,jp_tra) = ztrb(ji,jj,jk)
1108               psinkflx(ji,jj,jk)   = 2. * psinkflx(ji,jj,jk)
1109            END DO
1110         END DO
1111      END DO
1112!$OMP END DO NOWAIT
1113!$OMP END PARALLEL
1114      !
1115      CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
1116      !
1117      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink2')
1118      !
1119   END SUBROUTINE p4z_sink2
1120
1121
1122   INTEGER FUNCTION p4z_sink_alloc()
1123      !!----------------------------------------------------------------------
1124      !!                     ***  ROUTINE p4z_sink_alloc  ***
1125      !!----------------------------------------------------------------------
1126      ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4  (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) ,     &
1127         &      sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                      ,     &               
1128         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                      ,     &               
1129#if defined key_kriest
1130         &      xnumm(jpk)                                                        ,     &               
1131#else
1132         &      sinkfer2(jpi,jpj,jpk)                                             ,     &               
1133#endif
1134         &      sinkfer(jpi,jpj,jpk)                                              , STAT=p4z_sink_alloc )               
1135         !
1136      IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
1137      !
1138   END FUNCTION p4z_sink_alloc
1139   
1140#else
1141   !!======================================================================
1142   !!  Dummy module :                                   No PISCES bio-model
1143   !!======================================================================
1144CONTAINS
1145   SUBROUTINE p4z_sink                    ! Empty routine
1146   END SUBROUTINE p4z_sink
1147#endif 
1148
1149   !!======================================================================
1150END MODULE p4zsink
Note: See TracBrowser for help on using the repository browser.