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 @ 7508

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

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