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

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

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