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.
p5zsink.F90 in branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z – NEMO

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsink.F90 @ 6966

Last change on this file since 6966 was 6841, checked in by aumont, 8 years ago

Various bug fixes + explicit gamma function for lability
  • Property svn:executable set to *
File size: 33.9 KB
Line 
1MODULE p5zsink
2   !!======================================================================
3   !!                         ***  MODULE p5zsink  ***
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   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
11   !!----------------------------------------------------------------------
12#if defined key_pisces_quota
13   !!----------------------------------------------------------------------
14   !!   p5z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
15   !!   p5z_sink_init  :  Unitialisation of sinking speed parameters
16   !!   p5z_sink_alloc :  Allocate sinking speed variables
17   !!----------------------------------------------------------------------
18   USE oce_trc         !  shared variables between ocean and passive tracers
19   USE trc             !  passive tracers common variables
20   USE sms_pisces      !  PISCES Source Minus Sink variables
21   USE prtctl_trc      !  print control for debugging
22   USE iom             !  I/O manager
23   USE lib_mpp
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   p5z_sink         ! called in p5zbio.F90
29   PUBLIC   p5z_sink_init    ! called in trcsms_pisces.F90
30   PUBLIC   p5z_sink_alloc
31
32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio3   !: POC sinking speed
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio4   !: GOC sinking speed
34   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wscal    !: Calcite and BSi sinking speeds
35
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinking, sinking2  !: POC sinking fluxes
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingn, sinking2n  !: POC sinking fluxes
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkingp, sinking2p  !: POC sinking fluxes
39   !                                                          !  (different meanings depending on the parameterization)
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkcal, sinksil   !: CaCO3 and BSi sinking fluxes
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
42#if ! defined key_kriest
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
44#endif
45#if defined key_ligand
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfep      !: Fep sinking fluxes
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsfep
48#endif
49
50   INTEGER  :: ik100
51
52#if  defined key_kriest
53   REAL(wp) ::  xkr_sfact    !: Sinking factor
54   REAL(wp) ::  xkr_stick    !: Stickiness
55   REAL(wp) ::  xkr_nnano    !: Nbr of cell in nano size class
56   REAL(wp) ::  xkr_ndiat    !: Nbr of cell in diatoms size class
57   REAL(wp) ::  xkr_nmicro   !: Nbr of cell in microzoo size class
58   REAL(wp) ::  xkr_nmeso    !: Nbr of cell in mesozoo  size class
59   REAL(wp) ::  xkr_naggr    !: Nbr of cell in aggregates  size class
60
61   REAL(wp) ::  xkr_frac 
62
63   REAL(wp), PUBLIC ::  xkr_dnano       !: Size of particles in nano pool
64   REAL(wp), PUBLIC ::  xkr_ddiat       !: Size of particles in diatoms pool
65   REAL(wp), PUBLIC ::  xkr_dmicro      !: Size of particles in microzoo pool
66   REAL(wp), PUBLIC ::  xkr_dmeso       !: Size of particles in mesozoo pool
67   REAL(wp), PUBLIC ::  xkr_daggr       !: Size of particles in aggregates pool
68   REAL(wp), PUBLIC ::  xkr_wsbio_min   !: min vertical particle speed
69   REAL(wp), PUBLIC ::  xkr_wsbio_max   !: max vertical particle speed
70
71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   xnumm   !:  maximum number of particles in aggregates
72#endif
73
74   !!* Substitution
75#  include "top_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
78   !! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83#if ! defined key_kriest
84   !!----------------------------------------------------------------------
85   !!   'standard sinking parameterisation'                  ???
86   !!----------------------------------------------------------------------
87
88   SUBROUTINE p5z_sink ( kt, knt )
89      !!---------------------------------------------------------------------
90      !!                     ***  ROUTINE p5z_sink  ***
91      !!
92      !! ** Purpose :   Compute vertical flux of particulate matter due to
93      !!                gravitational sinking
94      !!
95      !! ** Method  : - ???
96      !!---------------------------------------------------------------------
97      INTEGER, INTENT(in) :: kt, knt
98      INTEGER  ::   ji, jj, jk, jit
99      INTEGER  ::   iiter1, iiter2
100      REAL(wp) ::   zfact, zwsmax, zmax
101      CHARACTER (len=25) :: charout
102      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
103      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
104      !!---------------------------------------------------------------------
105      !
106      IF( nn_timing == 1 )  CALL timing_start('p5z_sink')
107
108      ! Initialization of some global variables
109      ! ---------------------------------------
110      prodpoc(:,:,:) = 0.
111      conspoc(:,:,:) = 0.
112      prodgoc(:,:,:) = 0.
113      consgoc(:,:,:) = 0.
114      !
115      !    Sinking speeds of detritus is increased with depth as shown
116      !    by data and from the coagulation theory
117      !    -----------------------------------------------------------
118      DO jk = 1, jpkm1
119         DO jj = 1, jpj
120            DO ji = 1,jpi
121               zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
122               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / wsbio2scale
123               wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
124            END DO
125         END DO
126      END DO
127
128      ! limit the values of the sinking speeds to avoid numerical instabilities 
129      wsbio3(:,:,:) = wsbio
130#if defined key_ligand
131      wsfep (:,:,:) = wfep
132#endif
133      !
134      ! OA This is (I hope) a temporary solution for the problem that may
135      ! OA arise in specific situation where the CFL criterion is broken
136      ! OA for vertical sedimentation of particles. To avoid this, a time
137      ! OA splitting algorithm has been coded. A specific maximum
138      ! OA iteration number is provided and may be specified in the namelist
139      ! OA This is to avoid very large iteration number when explicit free
140      ! OA surface is used (for instance). When niter?max is set to 1,
141      ! OA this computation is skipped. The crude old threshold method is
142      ! OA then applied. This also happens when niter exceeds nitermax.
143      IF( MAX( niter1max, niter2max ) == 1 ) THEN
144        iiter1 = 1
145        iiter2 = 1
146      ELSE
147        iiter1 = 1
148        iiter2 = 1
149        DO jk = 1, jpkm1
150          DO jj = 1, jpj
151             DO ji = 1, jpi
152                IF( tmask(ji,jj,jk) == 1) THEN
153                   zwsmax =  0.5 * fse3t(ji,jj,jk) / xstep
154                   iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
155                   iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
156                ENDIF
157             END DO
158          END DO
159        END DO
160        IF( lk_mpp ) THEN
161           CALL mpp_max( iiter1 )
162           CALL mpp_max( iiter2 )
163        ENDIF
164        iiter1 = MIN( iiter1, niter1max )
165        iiter2 = MIN( iiter2, niter2max )
166      ENDIF
167
168      DO jk = 1,jpkm1
169         DO jj = 1, jpj
170            DO ji = 1, jpi
171               IF( tmask(ji,jj,jk) == 1 ) THEN
172                 zwsmax = 0.5 * fse3t(ji,jj,jk) / xstep
173                 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
174                 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
175#if defined key_ligand
176                 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
177#endif
178               ENDIF
179            END DO
180         END DO
181      END DO
182
183      wscal (:,:,:) = wsbio4(:,:,:)
184
185      !  Initializa to zero all the sinking arrays
186      !   -----------------------------------------
187      sinking (:,:,:) = 0.e0
188      sinking2(:,:,:) = 0.e0
189      sinkingn (:,:,:) = 0.e0
190      sinking2n(:,:,:) = 0.e0
191      sinkingp (:,:,:) = 0.e0
192      sinking2p(:,:,:) = 0.e0
193      sinkcal (:,:,:) = 0.e0
194      sinkfer (:,:,:) = 0.e0
195      sinksil (:,:,:) = 0.e0
196      sinkfer2(:,:,:) = 0.e0
197#if defined key_ligand
198      sinkfep(:,:,:) = 0.e0
199#endif
200
201      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
202      !   -----------------------------------------------------
203      DO jit = 1, iiter1
204        CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
205        CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 )
206        CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 )
207        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
208#if defined key_ligand
209        CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 )
210#endif
211      END DO
212
213      DO jit = 1, iiter2
214        CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
215        CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 )
216        CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 )
217        CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
218        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
219        CALL p4z_sink2( wscal, sinkcal , jpcal, iiter2 )
220      END DO
221
222     ! Total carbon export per year
223     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
224        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
225     !
226     IF( lk_iomput ) THEN
227       IF( knt == nrdttrc ) THEN
228          CALL wrk_alloc( jpi, jpj,      zw2d )
229          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
230          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
231          !
232          IF( iom_use( "EPC100" ) )  THEN
233              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
234              CALL iom_put( "EPC100"  , zw2d )
235          ENDIF
236          IF( iom_use( "EPFE100" ) )  THEN
237              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
238              CALL iom_put( "EPFE100"  , zw2d )
239          ENDIF
240          IF( iom_use( "EPCAL100" ) )  THEN
241              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
242              CALL iom_put( "EPCAL100"  , zw2d )
243          ENDIF
244          IF( iom_use( "EPSI100" ) )  THEN
245              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
246              CALL iom_put( "EPSI100"  , zw2d )
247          ENDIF
248          IF( iom_use( "EXPC" ) )  THEN
249              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
250              CALL iom_put( "EXPC"  , zw3d )
251          ENDIF
252          IF( iom_use( "EXPFE" ) )  THEN
253              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron
254              CALL iom_put( "EXPFE"  , zw3d )
255          ENDIF
256          IF( iom_use( "EXPCAL" ) )  THEN
257              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
258              CALL iom_put( "EXPCAL"  , zw3d )
259          ENDIF
260          IF( iom_use( "EXPSI" ) )  THEN
261              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
262              CALL iom_put( "EXPSI"  , zw3d )
263          ENDIF
264          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
265          !
266          CALL wrk_dealloc( jpi, jpj,      zw2d )
267          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
268        ENDIF
269      ELSE
270         IF( ln_diatrc ) THEN
271            zfact = 1.e3 * rfact2r
272            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
273            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
274            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
275            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)
276            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
277            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
278         ENDIF
279      ENDIF
280      !
281      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
282         WRITE(charout, FMT="('sink')")
283         CALL prt_ctl_trc_info(charout)
284         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
285      ENDIF
286      !
287      IF( nn_timing == 1 )  CALL timing_stop('p5z_sink')
288      !
289   END SUBROUTINE p5z_sink
290
291   SUBROUTINE p5z_sink_init
292      !!----------------------------------------------------------------------
293      !!                  ***  ROUTINE p5z_sink_init  ***
294      !!----------------------------------------------------------------------
295
296      INTEGER :: jk
297
298      ik100 = 10        !  last level where depth less than 100 m
299      DO jk = jpkm1, 1, -1
300         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
301      END DO
302      IF (lwp) WRITE(numout,*)
303      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
304      IF (lwp) WRITE(numout,*)
305      !
306      t_oce_co2_exp = 0._wp
307      !
308   END SUBROUTINE p5z_sink_init
309
310#else
311   !!----------------------------------------------------------------------
312   !!   'Kriest sinking parameterisation'        key_kriest          ???
313   !!----------------------------------------------------------------------
314
315   SUBROUTINE p5z_sink ( kt, knt )
316      !!---------------------------------------------------------------------
317      !!                ***  ROUTINE p5z_sink  ***
318      !!
319      !! ** Purpose :   Compute vertical flux of particulate matter due to
320      !!              gravitational sinking - Kriest parameterization
321      !!
322      !! ** Method  : - ???
323      !!---------------------------------------------------------------------
324      !
325      INTEGER, INTENT(in) :: kt, knt
326      !
327      INTEGER  :: ji, jj, jk, jit, niter1, niter2
328      REAL(wp) :: znum , zeps, zfm, zgm, zsm, zfactn, zfactp
329      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
330      REAL(wp) :: zval1, zval2, zval3, zval4
331      CHARACTER (len=25) :: charout
332      REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d 
333      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
334      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
335      !!---------------------------------------------------------------------
336      !
337      IF( nn_timing == 1 )  CALL timing_start('p5z_sink')
338      !
339      CALL wrk_alloc( jpi, jpj, jpk, znum3d )
340      !
341      !     Initialisation of variables used to compute Sinking Speed
342      !     ---------------------------------------------------------
343
344      znum3d(:,:,:) = 0.e0
345      zval1 = 1. + xkr_zeta
346      zval2 = 1. + xkr_zeta + xkr_eta
347      zval3 = 1. + xkr_eta
348
349      !     Computation of the vertical sinking speed : Kriest et Evans, 2000
350      !     -----------------------------------------------------------------
351
352      DO jk = 1, jpkm1
353         DO jj = 1, jpj
354            DO ji = 1, jpi
355               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
356                  znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
357                  ! -------------- To avoid sinking speed over 50 m/day -------
358                  znum  = MIN( xnumm(jk), znum )
359                  znum  = MAX( 1.1      , znum )
360                  znum3d(ji,jj,jk) = znum
361                  !------------------------------------------------------------
362                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
363                  zfm   = xkr_frac**( 1. - zeps )
364                  zgm   = xkr_frac**( zval1 - zeps )
365                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
366                  zdiv1 = zeps - zval3
367                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
368                     &             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
369                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
370                     &             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
371                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
372               ENDIF
373            END DO
374         END DO
375      END DO
376
377      wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )
378#if defined key_ligand
379      wsfep (:,:,:) = wfep
380#endif
381
382      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
383      !   -----------------------------------------
384
385      sinking (:,:,:) = 0.e0
386      sinkingn(:,:,:) = 0.e0
387      sinkingp(:,:,:) = 0.e0
388      sinking2(:,:,:) = 0.e0
389      sinkcal (:,:,:) = 0.e0
390      sinkfer (:,:,:) = 0.e0
391      sinksil (:,:,:) = 0.e0
392#if defined key_ligand
393      sinkfep(:,:,:) = 0.e0
394#endif
395
396     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
397     !   -----------------------------------------------------
398
399      niter1 = niter1max
400      niter2 = niter2max
401
402      DO jit = 1, niter1
403        CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )
404        CALL p4z_sink2( wsbio3, sinkingn, jppon, niter1 )
405        CALL p4z_sink2( wsbio3, sinkingp, jppop, niter1 )
406        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )
407        CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
408        CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
409#if defined key_ligand
410        CALL p4z_sink2( wsfep , sinkfep , jpfep, niter1 )
411#endif
412      END DO
413
414      DO jit = 1, niter2
415        CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )
416      END DO
417
418     ! Total carbon export per year
419     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
420        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
421     !
422     IF( lk_iomput ) THEN
423       IF( knt == nrdttrc ) THEN
424          CALL wrk_alloc( jpi, jpj,      zw2d )
425          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
426          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
427          !
428          IF( iom_use( "EPC100" ) )  THEN
429              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
430              CALL iom_put( "EPC100"  , zw2d )
431          ENDIF
432          IF( iom_use( "EPFE100" ) )  THEN
433              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
434              CALL iom_put( "EPFE100"  , zw2d )
435          ENDIF
436          IF( iom_use( "EPCAL100" ) )  THEN
437              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
438              CALL iom_put( "EPCAL100"  , zw2d )
439          ENDIF
440          IF( iom_use( "EPSI100" ) )  THEN
441              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
442              CALL iom_put( "EPSI100"  , zw2d )
443          ENDIF
444          IF( iom_use( "EXPC" ) )  THEN
445              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
446              CALL iom_put( "EXPC"  , zw3d )
447          ENDIF
448          IF( iom_use( "EXPFE" ) )  THEN
449              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron
450              CALL iom_put( "EXPFE"  , zw3d )
451          ENDIF
452          IF( iom_use( "EXPCAL" ) )  THEN
453              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
454              CALL iom_put( "EXPCAL"  , zw3d )
455          ENDIF
456          IF( iom_use( "EXPSI" ) )  THEN
457              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
458              CALL iom_put( "EXPSI"  , zw3d )
459          ENDIF
460          IF( iom_use( "XNUM" ) )  THEN
461              zw3d(:,:,:) =  znum3d(:,:,:) * tmask(:,:,:) !  Number of particles on aggregats
462              CALL iom_put( "XNUM"  , zw3d )
463          ENDIF
464          IF( iom_use( "WSC" ) )  THEN
465              zw3d(:,:,:) = wsbio3(:,:,:) * tmask(:,:,:) ! Sinking speed of carbon particles
466              CALL iom_put( "WSC"  , zw3d )
467          ENDIF
468          IF( iom_use( "WSN" ) )  THEN
469              zw3d(:,:,:) = wsbio4(:,:,:) * tmask(:,:,:) ! Sinking speed of particles number
470              CALL iom_put( "WSN"  , zw3d )
471          ENDIF
472          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
473          !
474          CALL wrk_dealloc( jpi, jpj,      zw2d )
475          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
476        ENDIF
477     ELSE
478        IF( ln_diatrc ) THEN
479            zfact = 1.e3 * rfact2r
480            trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik100) * zfact * tmask(:,:,1)
481            trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik100) * zfact * tmask(:,:,1)
482            trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik100) * zfact * tmask(:,:,1)
483            trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik100) * zfact * tmask(:,:,1)
484            trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik100) * zfact * tmask(:,:,1)
485            trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik100) * zfact * tmask(:,:,1)
486        ENDIF
487     ENDIF
488     !
489     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
490         WRITE(charout, FMT="('sink')")
491         CALL prt_ctl_trc_info(charout)
492         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
493     ENDIF
494     !
495     CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
496     !
497     IF( nn_timing == 1 )  CALL timing_stop('p5z_sink')
498     !
499   END SUBROUTINE p5z_sink
500
501
502   SUBROUTINE p5z_sink_init
503      !!----------------------------------------------------------------------
504      !!                  ***  ROUTINE p5z_sink_init  ***
505      !!
506      !! ** Purpose :   Initialization of sinking parameters
507      !!                Kriest parameterization only
508      !!
509      !! ** Method  :   Read the nampiskrs namelist and check the parameters
510      !!      called at the first timestep
511      !!
512      !! ** input   :   Namelist nampiskrs
513      !!----------------------------------------------------------------------
514      INTEGER  ::   jk, jn, kiter
515      INTEGER  ::   ios                 ! Local integer output status for namelist read
516      REAL(wp) ::   znum, zdiv
517      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
518      REAL(wp) ::   zmin, zmax, zl, zr, xacc
519      !
520      NAMELIST/nampiskrs/ xkr_sfact, xkr_stick ,  &
521         &                xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr
522      !!----------------------------------------------------------------------
523      !
524      IF( nn_timing == 1 )  CALL timing_start('p5z_sink_init')
525      !
526
527      REWIND( numnatp_ref )              ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest
528      READ  ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)
529901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )
530
531      REWIND( numnatp_cfg )              ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest
532      READ  ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )
533902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )
534      IF(lwm) WRITE ( numonp, nampiskrs )
535
536      IF(lwp) THEN
537         WRITE(numout,*)
538         WRITE(numout,*) ' Namelist : nampiskrs'
539         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
540         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
541         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
542         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
543         WRITE(numout,*) '    Nbr of cell in microzoo size class       xkr_nmicro   = ', xkr_nmicro
544         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
545         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
546      ENDIF
547
548
549      ! max and min vertical particle speed
550      xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
551      xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
552      IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
553
554      !
555      !    effect of the sizes of the different living pools on particle numbers
556      !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
557      !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
558      !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
559      !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
560      !    doc aggregates = 1um
561      ! ----------------------------------------------------------
562
563      xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
564      xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
565      xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )
566      xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
567      xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
568
569      !!---------------------------------------------------------------------
570      !!    'key_kriest'                                                  ???
571      !!---------------------------------------------------------------------
572      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
573      !  Search of the maximum number of particles in aggregates for each k-level.
574      !  Bissection Method
575      !--------------------------------------------------------------------
576      IF (lwp) THEN
577        WRITE(numout,*)
578        WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
579      ENDIF
580
581      xacc     =  0.001_wp
582      kiter    = 50
583      zmin     =  1.10_wp
584      zmax     = xkr_mass_max / xkr_mass_min
585      xkr_frac = zmax
586
587      DO jk = 1,jpk
588         zl = zmin
589         zr = zmax
590         wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2
591         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
592         znum = zl - 1.
593         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
594            & - ( xkr_wsbio_max * xkr_eta * znum * &
595            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
596            & - wmax
597
598         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
599         znum = zr - 1.
600         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
601            & - ( xkr_wsbio_max * xkr_eta * znum * &
602            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
603            & - wmax
604iflag:   DO jn = 1, kiter
605            IF    ( zwl == 0._wp ) THEN   ;   znummax = zl
606            ELSEIF( zwr == 0._wp ) THEN   ;   znummax = zr
607            ELSE
608               znummax = ( zr + zl ) / 2.
609               zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
610               znum = znummax - 1.
611               zws =  xkr_wsbio_min * xkr_zeta / zdiv &
612                  & - ( xkr_wsbio_max * xkr_eta * znum * &
613                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
614                  & - wmax
615               IF( zws * zwl < 0. ) THEN   ;   zr = znummax
616               ELSE                        ;   zl = znummax
617               ENDIF
618               zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
619               znum = zl - 1.
620               zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
621                  & - ( xkr_wsbio_max * xkr_eta * znum * &
622                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
623                  & - wmax
624
625               zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
626               znum = zr - 1.
627               zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
628                  & - ( xkr_wsbio_max * xkr_eta * znum * &
629                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
630                  & - wmax
631               !
632               IF ( ABS ( zws )  <= xacc ) EXIT iflag
633               !
634            ENDIF
635            !
636         END DO iflag
637
638         xnumm(jk) = znummax
639         IF (lwp) WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
640         !
641      END DO
642      !
643      ik100 = 10        !  last level where depth less than 100 m
644      DO jk = jpkm1, 1, -1
645         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
646      END DO
647      IF (lwp) WRITE(numout,*)
648      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
649      IF (lwp) WRITE(numout,*)
650      !
651      t_oce_co2_exp = 0._wp
652      !
653      IF( nn_timing == 1 )  CALL timing_stop('p5z_sink_init')
654      !
655  END SUBROUTINE p5z_sink_init
656
657#endif
658
659   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
660      !!---------------------------------------------------------------------
661      !!                     ***  ROUTINE p4z_sink2  ***
662      !!
663      !! ** Purpose :   Compute the sedimentation terms for the various sinking
664      !!     particles. The scheme used to compute the trends is based
665      !!     on MUSCL.
666      !!
667      !! ** Method  : - this ROUTINE compute not exactly the advection but the
668      !!      transport term, i.e.  div(u*tra).
669      !!---------------------------------------------------------------------
670      !
671      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
672      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting
673      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
674      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
675      !!
676      INTEGER  ::   ji, jj, jk, jn
677      REAL(wp) ::   zigma,zew,zign, zflx, zstep
678      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb 
679      !!---------------------------------------------------------------------
680      !
681      IF( nn_timing == 1 )  CALL timing_start('p4z_sink2')
682      !
683      ! Allocate temporary workspace
684      CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
685
686      zstep = rfact2 / FLOAT( kiter ) / 2.
687
688      ztraz(:,:,:) = 0.e0
689      zakz (:,:,:) = 0.e0
690      ztrb (:,:,:) = trb(:,:,:,jp_tra)
691
692      DO jk = 1, jpkm1
693         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
694      END DO
695      zwsink2(:,:,1) = 0.e0
696      IF( lk_degrad ) THEN
697         zwsink2(:,:,:) = zwsink2(:,:,:) * facvol(:,:,:)
698      ENDIF
699
700
701      ! Vertical advective flux
702      DO jn = 1, 2
703         !  first guess of the slopes interior values
704         DO jk = 2, jpkm1
705            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
706         END DO
707         ztraz(:,:,1  ) = 0.0
708         ztraz(:,:,jpk) = 0.0
709
710         ! slopes
711         DO jk = 2, jpkm1
712            DO jj = 1,jpj
713               DO ji = 1, jpi
714                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
715                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
716               END DO
717            END DO
718         END DO
719         
720         ! Slopes limitation
721         DO jk = 2, jpkm1
722            DO jj = 1, jpj
723               DO ji = 1, jpi
724                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
725                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
726               END DO
727            END DO
728         END DO
729         
730         ! vertical advective flux
731         DO jk = 1, jpkm1
732            DO jj = 1, jpj     
733               DO ji = 1, jpi   
734                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
735                  zew   = zwsink2(ji,jj,jk+1)
736                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
737               END DO
738            END DO
739         END DO
740         !
741         ! Boundary conditions
742         psinkflx(:,:,1  ) = 0.e0
743         psinkflx(:,:,jpk) = 0.e0
744         
745         DO jk=1,jpkm1
746            DO jj = 1,jpj
747               DO ji = 1, jpi
748                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
749                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
750               END DO
751            END DO
752         END DO
753
754      ENDDO
755
756      DO jk = 1,jpkm1
757         DO jj = 1,jpj
758            DO ji = 1, jpi
759               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
760               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
761            END DO
762         END DO
763      END DO
764
765      trb(:,:,:,jp_tra) = ztrb(:,:,:)
766      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
767      !
768      CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
769      !
770      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink2')
771      !
772   END SUBROUTINE p4z_sink2
773
774   INTEGER FUNCTION p5z_sink_alloc()
775      !!----------------------------------------------------------------------
776      !!                     ***  ROUTINE p5z_sink_alloc  ***
777      !!----------------------------------------------------------------------
778      ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4  (jpi,jpj,jpk)   , wscal(jpi,jpj,jpk) ,     &
779         &      sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)   ,     &               
780         &      sinkingn(jpi,jpj,jpk) , sinking2n(jpi,jpj,jpk) ,     &
781         &      sinkingp(jpi,jpj,jpk) , sinking2p(jpi,jpj,jpk) ,     &
782
783         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)   ,     &   
784#if defined key_kriest
785         &      xnumm(jpk)                                                        ,     &               
786#else
787         &      sinkfer2(jpi,jpj,jpk)                                             ,     &               
788#endif
789#if defined key_ligand
790         &      wsfep(jpi,jpj,jpk)  , sinkfep(jpi,jpj,jpk)     ,     &
791#endif
792         &      sinkfer(jpi,jpj,jpk)     , STAT=p5z_sink_alloc )               
793         !
794      IF( p5z_sink_alloc /= 0 ) CALL ctl_warn('p5z_sink_alloc : failed to allocate arrays.')
795      !
796   END FUNCTION p5z_sink_alloc
797   
798#else
799   !!======================================================================
800   !!  Dummy module :                                   No PISCES bio-model
801   !!======================================================================
802CONTAINS
803   SUBROUTINE p5z_sink                    ! Empty routine
804   END SUBROUTINE p5z_sink
805#endif 
806
807   !!======================================================================
808END MODULE  p5zsink
Note: See TracBrowser for help on using the repository browser.