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/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
File size: 40.5 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   !! * Substitutions
68#  include "domzgr_substitute.h90"
69   !!----------------------------------------------------------------------
70   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
71   !! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
72   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74CONTAINS
75
76#if ! defined key_kriest
77   !!----------------------------------------------------------------------
78   !!   'standard sinking parameterisation'                  ???
79   !!----------------------------------------------------------------------
80
81   SUBROUTINE p4z_sink ( kt, knt )
82      !!---------------------------------------------------------------------
83      !!                     ***  ROUTINE p4z_sink  ***
84      !!
85      !! ** Purpose :   Compute vertical flux of particulate matter due to
86      !!                gravitational sinking
87      !!
88      !! ** Method  : - ???
89      !!---------------------------------------------------------------------
90      INTEGER, INTENT(in) :: kt, knt
91      INTEGER  ::   ji, jj, jk, jit
92      INTEGER  ::   iiter1, iiter2
93      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
94      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
95      REAL(wp) ::   zfact, zwsmax, zmax, zstep
96      CHARACTER (len=25) :: charout
97      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
98      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
99      !!---------------------------------------------------------------------
100      !
101      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
102      !
103      !    Sinking speeds of detritus is increased with depth as shown
104      !    by data and from the coagulation theory
105      !    -----------------------------------------------------------
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., fsdepw(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      wsbio3(:,:,:) = wsbio
118      wscal (:,:,:) = wsbio4(:,:,:)
119      !
120      ! OA This is (I hope) a temporary solution for the problem that may
121      ! OA arise in specific situation where the CFL criterion is broken
122      ! OA for vertical sedimentation of particles. To avoid this, a time
123      ! OA splitting algorithm has been coded. A specific maximum
124      ! OA iteration number is provided and may be specified in the namelist
125      ! OA This is to avoid very large iteration number when explicit free
126      ! OA surface is used (for instance). When niter?max is set to 1,
127      ! OA this computation is skipped. The crude old threshold method is
128      ! OA then applied. This also happens when niter exceeds nitermax.
129      IF( MAX( niter1max, niter2max ) == 1 ) THEN
130        iiter1 = 1
131        iiter2 = 1
132      ELSE
133        iiter1 = 1
134        iiter2 = 1
135        DO jk = 1, jpkm1
136          DO jj = 1, jpj
137             DO ji = 1, jpi
138                IF( tmask(ji,jj,jk) == 1) THEN
139                   zwsmax =  0.5 * fse3t(ji,jj,jk) / xstep
140                   iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
141                   iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
142                ENDIF
143             END DO
144          END DO
145        END DO
146        IF( lk_mpp ) THEN
147           CALL mpp_max( iiter1 )
148           CALL mpp_max( iiter2 )
149        ENDIF
150        iiter1 = MIN( iiter1, niter1max )
151        iiter2 = MIN( iiter2, niter2max )
152      ENDIF
153
154      DO jk = 1,jpkm1
155         DO jj = 1, jpj
156            DO ji = 1, jpi
157               IF( tmask(ji,jj,jk) == 1 ) THEN
158                 zwsmax = 0.5 * fse3t(ji,jj,jk) / xstep
159                 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
160                 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
161               ENDIF
162            END DO
163         END DO
164      END DO
165
166      !  Initializa to zero all the sinking arrays
167      !   -----------------------------------------
168      sinking (:,:,:) = 0.e0
169      sinking2(:,:,:) = 0.e0
170      sinkcal (:,:,:) = 0.e0
171      sinkfer (:,:,:) = 0.e0
172      sinksil (:,:,:) = 0.e0
173      sinkfer2(:,:,:) = 0.e0
174
175      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
176      !   -----------------------------------------------------
177      DO jit = 1, iiter1
178        CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
179        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
180      END DO
181
182      DO jit = 1, iiter2
183        CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
184        CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
185        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
186        CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
187      END DO
188
189      !  Exchange between organic matter compartments due to coagulation/disaggregation
190      !  ---------------------------------------------------
191      DO jk = 1, jpkm1
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194               !
195               zstep = xstep 
196# if defined key_degrad
197               zstep = zstep * facvol(ji,jj,jk)
198# endif
199               zfact = zstep * xdiss(ji,jj,jk)
200               !  Part I : Coagulation dependent on turbulence
201               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
202               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
203
204               ! Part II : Differential settling
205
206               !  Aggregation of small into large particles
207               zagg3 =  47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
208               zagg4 =  3.3  * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
209
210               zagg   = zagg1 + zagg2 + zagg3 + zagg4
211               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
212
213               ! Aggregation of DOC to POC :
214               ! 1st term is shear aggregation of DOC-DOC
215               ! 2nd term is shear aggregation of DOC-POC
216               ! 3rd term is differential settling of DOC-POC
217               zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
218               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
219               ! transfer of DOC to GOC :
220               ! 1st term is shear aggregation
221               ! 2nd term is differential settling
222               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
223               ! tranfer of DOC to POC due to brownian motion
224               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
225
226               !  Update the trends
227               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
228               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
229               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
230               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
231               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
232               !
233            END DO
234         END DO
235      END DO
236
237
238     ! Total carbon export per year
239     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
240        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
241     !
242     IF( lk_iomput ) THEN
243       IF( knt == nrdttrc ) THEN
244          CALL wrk_alloc( jpi, jpj,      zw2d )
245          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
246          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
247          !
248          IF( iom_use( "EPC100" ) )  THEN
249              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
250              CALL iom_put( "EPC100"  , zw2d )
251          ENDIF
252          IF( iom_use( "EPFE100" ) )  THEN
253              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
254              CALL iom_put( "EPFE100"  , zw2d )
255          ENDIF
256          IF( iom_use( "EPCAL100" ) )  THEN
257              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
258              CALL iom_put( "EPCAL100"  , zw2d )
259          ENDIF
260          IF( iom_use( "EPSI100" ) )  THEN
261              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
262              CALL iom_put( "EPSI100"  , zw2d )
263          ENDIF
264          IF( iom_use( "EXPC" ) )  THEN
265              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
266              CALL iom_put( "EXPC"  , zw3d )
267          ENDIF
268          IF( iom_use( "EXPFE" ) )  THEN
269              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron
270              CALL iom_put( "EXPFE"  , zw3d )
271          ENDIF
272          IF( iom_use( "EXPCAL" ) )  THEN
273              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
274              CALL iom_put( "EXPCAL"  , zw3d )
275          ENDIF
276          IF( iom_use( "EXPSI" ) )  THEN
277              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
278              CALL iom_put( "EXPSI"  , zw3d )
279          ENDIF
280          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
281          !
282          CALL wrk_dealloc( jpi, jpj,      zw2d )
283          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
284        ENDIF
285      ELSE
286         IF( ln_diatrc ) THEN
287            zfact = 1.e3 * rfact2r
288            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
289            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
290            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
291            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)
292            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
293            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
294         ENDIF
295      ENDIF
296      !
297      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
298         WRITE(charout, FMT="('sink')")
299         CALL prt_ctl_trc_info(charout)
300         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
301      ENDIF
302      !
303      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
304      !
305   END SUBROUTINE p4z_sink
306
307   SUBROUTINE p4z_sink_init
308      !!----------------------------------------------------------------------
309      !!                  ***  ROUTINE p4z_sink_init  ***
310      !!----------------------------------------------------------------------
311      INTEGER :: jk
312
313      ik100 = 10        !  last level where depth less than 100 m
314      DO jk = jpkm1, 1, -1
315         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
316      END DO
317      IF (lwp) WRITE(numout,*)
318      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
319      IF (lwp) WRITE(numout,*)
320      !
321      t_oce_co2_exp = 0._wp
322      !
323   END SUBROUTINE p4z_sink_init
324
325#else
326   !!----------------------------------------------------------------------
327   !!   'Kriest sinking parameterisation'        key_kriest          ???
328   !!----------------------------------------------------------------------
329
330   SUBROUTINE p4z_sink ( kt, knt )
331      !!---------------------------------------------------------------------
332      !!                ***  ROUTINE p4z_sink  ***
333      !!
334      !! ** Purpose :   Compute vertical flux of particulate matter due to
335      !!              gravitational sinking - Kriest parameterization
336      !!
337      !! ** Method  : - ???
338      !!---------------------------------------------------------------------
339      !
340      INTEGER, INTENT(in) :: kt, knt
341      !
342      INTEGER  :: ji, jj, jk, jit, niter1, niter2
343      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh
344      REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc
345      REAL(wp) :: znum , zeps, zfm, zgm, zsm
346      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
347      REAL(wp) :: zval1, zval2, zval3, zval4
348      REAL(wp) :: zfact
349      INTEGER  :: ik1
350      CHARACTER (len=25) :: charout
351      REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d 
352      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
353      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
354      !!---------------------------------------------------------------------
355      !
356      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
357      !
358      CALL wrk_alloc( jpi, jpj, jpk, znum3d )
359      !
360      !     Initialisation of variables used to compute Sinking Speed
361      !     ---------------------------------------------------------
362
363      znum3d(:,:,:) = 0.e0
364      zval1 = 1. + xkr_zeta
365      zval2 = 1. + xkr_zeta + xkr_eta
366      zval3 = 1. + xkr_eta
367
368      !     Computation of the vertical sinking speed : Kriest et Evans, 2000
369      !     -----------------------------------------------------------------
370
371      DO jk = 1, jpkm1
372         DO jj = 1, jpj
373            DO ji = 1, jpi
374               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
375                  znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
376                  ! -------------- To avoid sinking speed over 50 m/day -------
377                  znum  = MIN( xnumm(jk), znum )
378                  znum  = MAX( 1.1      , znum )
379                  znum3d(ji,jj,jk) = znum
380                  !------------------------------------------------------------
381                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
382                  zfm   = xkr_frac**( 1. - zeps )
383                  zgm   = xkr_frac**( zval1 - zeps )
384                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
385                  zdiv1 = zeps - zval3
386                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
387                     &             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
388                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
389                     &             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
390                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
391               ENDIF
392            END DO
393         END DO
394      END DO
395
396      wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )
397
398      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
399      !   -----------------------------------------
400
401      sinking (:,:,:) = 0.e0
402      sinking2(:,:,:) = 0.e0
403      sinkcal (:,:,:) = 0.e0
404      sinkfer (:,:,:) = 0.e0
405      sinksil (:,:,:) = 0.e0
406
407     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
408     !   -----------------------------------------------------
409
410      niter1 = niter1max
411      niter2 = niter2max
412
413      DO jit = 1, niter1
414        CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )
415        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )
416        CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
417        CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
418      END DO
419
420      DO jit = 1, niter2
421        CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )
422      END DO
423
424     !  Exchange between organic matter compartments due to coagulation/disaggregation
425     !  ---------------------------------------------------
426
427      zval1 = 1. + xkr_zeta
428      zval2 = 1. + xkr_eta
429      zval3 = 3. + xkr_eta
430      zval4 = 4. + xkr_eta
431
432      DO jk = 1,jpkm1
433         DO jj = 1,jpj
434            DO ji = 1,jpi
435               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
436
437                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp
438                  !-------------- To avoid sinking speed over 50 m/day -------
439                  znum  = min(xnumm(jk),znum)
440                  znum  = MAX( 1.1,znum)
441                  !------------------------------------------------------------
442                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
443                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
444                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
445                  zdiv2 = zeps - 2.
446                  zdiv3 = zeps - 3.
447                  zdiv4 = zeps - zval2
448                  zdiv5 = 2.* zeps - zval4
449                  zfm   = xkr_frac**( 1.- zeps )
450                  zsm   = xkr_frac**xkr_eta
451
452                  !    Part I : Coagulation dependant on turbulence
453                  !    ----------------------------------------------
454
455                  zagg1 =  0.163 * trb(ji,jj,jk,jpnum)**2               &
456                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
457                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
458                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
459                     &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
460                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm*                       &
461                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
462                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
463                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
464                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
465                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
466                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))   
467
468                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3 
469                 
470                 !    Aggregation of small into large particles
471                 !    Part II : Differential settling
472                 !    ----------------------------------------------
473
474                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       &
475                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
476                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
477                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
478                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
479                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
480
481                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         &
482                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
483                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
484                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
485                     &                 /zdiv) 
486
487                  !
488                  !     Fractionnation by swimming organisms
489                  !     ------------------------------------
490
491                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  &
492                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
493                    &      * 10000.*xstep
494
495                  !     Aggregation of DOC to small particles
496                  !     --------------------------------------
497
498                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   &
499                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)
500                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  &
501                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)
502
503# if defined key_degrad
504                   zagg1   = zagg1   * facvol(ji,jj,jk)                 
505                   zagg2   = zagg2   * facvol(ji,jj,jk)                 
506                   zagg3   = zagg3   * facvol(ji,jj,jk)                 
507                   zagg4   = zagg4   * facvol(ji,jj,jk)                 
508                   zagg5   = zagg5   * facvol(ji,jj,jk)                 
509                   zaggdoc = zaggdoc * facvol(ji,jj,jk)                 
510                   zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)
511# endif
512                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
513                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
514                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
515                  !
516                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
517                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
518                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
519                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1
520
521               ENDIF
522            END DO
523         END DO
524      END DO
525
526     ! Total primary production per year
527     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,ik100) * e1e2t(:,:) * tmask(:,:,1) )
528     !
529     IF( lk_iomput ) THEN
530        IF( knt == nrdttrc ) THEN
531          CALL wrk_alloc( jpi, jpj,      zw2d )
532          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
533          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
534          !
535          IF( iom_use( "EPC100" ) )  THEN
536              zw2d(:,:) = sinking(:,:,ik100) * zfact * tmask(:,:,1) ! Export of carbon at 100m
537              CALL iom_put( "EPC100"  , zw2d )
538          ENDIF
539          IF( iom_use( "EPN100" ) )  THEN
540              zw2d(:,:) = sinking2(:,:,ik100) * zfact * tmask(:,:,1) ! Export of number of aggregates ?
541              CALL iom_put( "EPN100"  , zw2d )
542          ENDIF
543          IF( iom_use( "EPCAL100" ) )  THEN
544              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
545              CALL iom_put( "EPCAL100"  , zw2d )
546          ENDIF
547          IF( iom_use( "EPSI100" ) )  THEN
548              zw2d(:,:) = sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
549              CALL iom_put( "EPSI100"  , zw2d )
550          ENDIF
551          IF( iom_use( "EXPC" ) )  THEN
552              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
553              CALL iom_put( "EXPC"  , zw3d )
554          ENDIF
555          IF( iom_use( "EXPN" ) )  THEN
556              zw3d(:,:,:) = sinking(:,:,:) * zfact * tmask(:,:,:) ! Export of carbon in the water column
557              CALL iom_put( "EXPN"  , zw3d )
558          ENDIF
559          IF( iom_use( "EXPCAL" ) )  THEN
560              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
561              CALL iom_put( "EXPCAL"  , zw3d )
562          ENDIF
563          IF( iom_use( "EXPSI" ) )  THEN
564              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
565              CALL iom_put( "EXPSI"  , zw3d )
566          ENDIF
567          IF( iom_use( "XNUM" ) )  THEN
568              zw3d(:,:,:) =  znum3d(:,:,:) * tmask(:,:,:) !  Number of particles on aggregats
569              CALL iom_put( "XNUM"  , zw3d )
570          ENDIF
571          IF( iom_use( "WSC" ) )  THEN
572              zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles
573              CALL iom_put( "WSC"  , zw3d )
574          ENDIF
575          IF( iom_use( "WSN" ) )  THEN
576              zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number
577              CALL iom_put( "WSN"  , zw3d )
578          ENDIF
579          !
580          CALL wrk_dealloc( jpi, jpj,      zw2d )
581          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
582      ELSE
583         IF( ln_diatrc ) THEN
584            zfact = 1.e3 * rfact2r
585            trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik100)  * zfact * tmask(:,:,1)
586            trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik100)  * zfact * tmask(:,:,1)
587            trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik100)  * zfact * tmask(:,:,1)
588            trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik100)  * zfact * tmask(:,:,1)
589            trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik100)  * zfact * tmask(:,:,1)
590            trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zfact * tmask(:,:,:)
591            trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zfact * tmask(:,:,:)
592            trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zfact * tmask(:,:,:)
593            trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zfact * tmask(:,:,:)
594            trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)              * tmask(:,:,:)
595            trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)              * tmask(:,:,:)
596            trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)              * tmask(:,:,:)
597         ENDIF
598      ENDIF
599
600      !
601      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
602         WRITE(charout, FMT="('sink')")
603         CALL prt_ctl_trc_info(charout)
604         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
605      ENDIF
606      !
607      CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
608      !
609      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
610      !
611   END SUBROUTINE p4z_sink
612
613
614   SUBROUTINE p4z_sink_init
615      !!----------------------------------------------------------------------
616      !!                  ***  ROUTINE p4z_sink_init  ***
617      !!
618      !! ** Purpose :   Initialization of sinking parameters
619      !!                Kriest parameterization only
620      !!
621      !! ** Method  :   Read the nampiskrs namelist and check the parameters
622      !!      called at the first timestep
623      !!
624      !! ** input   :   Namelist nampiskrs
625      !!----------------------------------------------------------------------
626      INTEGER  ::   jk, jn, kiter
627      INTEGER  ::   ios                 ! Local integer output status for namelist read
628      REAL(wp) ::   znum, zdiv
629      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
630      REAL(wp) ::   zmin, zmax, zl, zr, xacc
631      !
632      NAMELIST/nampiskrs/ xkr_sfact, xkr_stick ,  &
633         &                xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr
634      !!----------------------------------------------------------------------
635      !
636      IF( nn_timing == 1 )  CALL timing_start('p4z_sink_init')
637      !
638
639      REWIND( numnatp_ref )              ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest
640      READ  ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)
641901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )
642
643      REWIND( numnatp_cfg )              ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest
644      READ  ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )
645902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )
646      IF(lwm) WRITE ( numonp, nampiskrs )
647
648      IF(lwp) THEN
649         WRITE(numout,*)
650         WRITE(numout,*) ' Namelist : nampiskrs'
651         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
652         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
653         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
654         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
655         WRITE(numout,*) '    Nbr of cell in microzoo size class       xkr_nmicro   = ', xkr_nmicro
656         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
657         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
658      ENDIF
659
660
661      ! max and min vertical particle speed
662      xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
663      xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
664      IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
665
666      !
667      !    effect of the sizes of the different living pools on particle numbers
668      !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
669      !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
670      !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
671      !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
672      !    doc aggregates = 1um
673      ! ----------------------------------------------------------
674
675      xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
676      xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
677      xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )
678      xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
679      xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
680
681      !!---------------------------------------------------------------------
682      !!    'key_kriest'                                                  ???
683      !!---------------------------------------------------------------------
684      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
685      !  Search of the maximum number of particles in aggregates for each k-level.
686      !  Bissection Method
687      !--------------------------------------------------------------------
688      IF (lwp) THEN
689        WRITE(numout,*)
690        WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
691      ENDIF
692
693      xacc     =  0.001_wp
694      kiter    = 50
695      zmin     =  1.10_wp
696      zmax     = xkr_mass_max / xkr_mass_min
697      xkr_frac = zmax
698
699      DO jk = 1,jpk
700         zl = zmin
701         zr = zmax
702         wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2
703         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
704         znum = zl - 1.
705         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
706            & - ( xkr_wsbio_max * xkr_eta * znum * &
707            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
708            & - wmax
709
710         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
711         znum = zr - 1.
712         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
713            & - ( xkr_wsbio_max * xkr_eta * znum * &
714            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
715            & - wmax
716iflag:   DO jn = 1, kiter
717            IF    ( zwl == 0._wp ) THEN   ;   znummax = zl
718            ELSEIF( zwr == 0._wp ) THEN   ;   znummax = zr
719            ELSE
720               znummax = ( zr + zl ) / 2.
721               zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
722               znum = znummax - 1.
723               zws =  xkr_wsbio_min * xkr_zeta / zdiv &
724                  & - ( xkr_wsbio_max * xkr_eta * znum * &
725                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
726                  & - wmax
727               IF( zws * zwl < 0. ) THEN   ;   zr = znummax
728               ELSE                        ;   zl = znummax
729               ENDIF
730               zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
731               znum = zl - 1.
732               zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
733                  & - ( xkr_wsbio_max * xkr_eta * znum * &
734                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
735                  & - wmax
736
737               zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
738               znum = zr - 1.
739               zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
740                  & - ( xkr_wsbio_max * xkr_eta * znum * &
741                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
742                  & - wmax
743               !
744               IF ( ABS ( zws )  <= xacc ) EXIT iflag
745               !
746            ENDIF
747            !
748         END DO iflag
749
750         xnumm(jk) = znummax
751         IF (lwp) WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
752         !
753      END DO
754      !
755      ik100 = 10        !  last level where depth less than 100 m
756      DO jk = jpkm1, 1, -1
757         IF( gdept_1d(jk) > 100. )  iksed = jk - 1
758      END DO
759      IF (lwp) WRITE(numout,*)
760      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
761      IF (lwp) WRITE(numout,*)
762      !
763      t_oce_co2_exp = 0._wp
764      !
765      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink_init')
766      !
767  END SUBROUTINE p4z_sink_init
768
769#endif
770
771   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
772      !!---------------------------------------------------------------------
773      !!                     ***  ROUTINE p4z_sink2  ***
774      !!
775      !! ** Purpose :   Compute the sedimentation terms for the various sinking
776      !!     particles. The scheme used to compute the trends is based
777      !!     on MUSCL.
778      !!
779      !! ** Method  : - this ROUTINE compute not exactly the advection but the
780      !!      transport term, i.e.  div(u*tra).
781      !!---------------------------------------------------------------------
782      !
783      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
784      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting
785      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
786      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
787      !!
788      INTEGER  ::   ji, jj, jk, jn
789      REAL(wp) ::   zigma,zew,zign, zflx, zstep
790      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb 
791      !!---------------------------------------------------------------------
792      !
793      IF( nn_timing == 1 )  CALL timing_start('p4z_sink2')
794      !
795      ! Allocate temporary workspace
796      CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
797
798      zstep = rfact2 / FLOAT( kiter ) / 2.
799
800      ztraz(:,:,:) = 0.e0
801      zakz (:,:,:) = 0.e0
802      ztrb (:,:,:) = trb(:,:,:,jp_tra)
803
804      DO jk = 1, jpkm1
805         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
806      END DO
807      zwsink2(:,:,1) = 0.e0
808      IF( lk_degrad ) THEN
809         zwsink2(:,:,:) = zwsink2(:,:,:) * facvol(:,:,:)
810      ENDIF
811
812
813      ! Vertical advective flux
814      DO jn = 1, 2
815         !  first guess of the slopes interior values
816         DO jk = 2, jpkm1
817            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
818         END DO
819         ztraz(:,:,1  ) = 0.0
820         ztraz(:,:,jpk) = 0.0
821
822         ! slopes
823         DO jk = 2, jpkm1
824            DO jj = 1,jpj
825               DO ji = 1, jpi
826                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
827                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
828               END DO
829            END DO
830         END DO
831         
832         ! Slopes limitation
833         DO jk = 2, jpkm1
834            DO jj = 1, jpj
835               DO ji = 1, jpi
836                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
837                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
838               END DO
839            END DO
840         END DO
841         
842         ! vertical advective flux
843         DO jk = 1, jpkm1
844            DO jj = 1, jpj     
845               DO ji = 1, jpi   
846                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
847                  zew   = zwsink2(ji,jj,jk+1)
848                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
849               END DO
850            END DO
851         END DO
852         !
853         ! Boundary conditions
854         psinkflx(:,:,1  ) = 0.e0
855         psinkflx(:,:,jpk) = 0.e0
856         
857         DO jk=1,jpkm1
858            DO jj = 1,jpj
859               DO ji = 1, jpi
860                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
861                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
862               END DO
863            END DO
864         END DO
865
866      ENDDO
867
868      DO jk = 1,jpkm1
869         DO jj = 1,jpj
870            DO ji = 1, jpi
871               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
872               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
873            END DO
874         END DO
875      END DO
876
877      trb(:,:,:,jp_tra) = ztrb(:,:,:)
878      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
879      !
880      CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
881      !
882      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink2')
883      !
884   END SUBROUTINE p4z_sink2
885
886
887   INTEGER FUNCTION p4z_sink_alloc()
888      !!----------------------------------------------------------------------
889      !!                     ***  ROUTINE p4z_sink_alloc  ***
890      !!----------------------------------------------------------------------
891      ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4  (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) ,     &
892         &      sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                      ,     &               
893         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                      ,     &               
894#if defined key_kriest
895         &      xnumm(jpk)                                                        ,     &               
896#else
897         &      sinkfer2(jpi,jpj,jpk)                                             ,     &               
898#endif
899         &      sinkfer(jpi,jpj,jpk)                                              , STAT=p4z_sink_alloc )               
900         !
901      IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
902      !
903   END FUNCTION p4z_sink_alloc
904   
905#else
906   !!======================================================================
907   !!  Dummy module :                                   No PISCES bio-model
908   !!======================================================================
909CONTAINS
910   SUBROUTINE p4z_sink                    ! Empty routine
911   END SUBROUTINE p4z_sink
912#endif 
913
914   !!======================================================================
915END MODULE p4zsink
Note: See TracBrowser for help on using the repository browser.