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

source: trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 4641

Last change on this file since 4641 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
File size: 36.4 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  :: iksed  = 10
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   !!* Substitution
68#  include "top_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, jnt )
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, jnt
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      REAL(wp) ::   zrfact2
97      INTEGER  ::   ik1
98      CHARACTER (len=25) :: charout
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 * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
202               zagg2 = 4452. * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
203
204               ! Part II : Differential settling
205
206               !  Aggregation of small into large particles
207               zagg3 =  47.1 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
208               zagg4 =  3.3  * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
209
210               zagg   = zagg1 + zagg2 + zagg3 + zagg4
211               zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(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 * trn(ji,jj,jk,jpdoc) + 102.4 * trn(ji,jj,jk,jppoc) ) * zfact       &
218               &            + 2.4 * zstep * trn(ji,jj,jk,jppoc) ) * 0.3 * trn(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 ) * trn(ji,jj,jk,jpgoc) * 0.3 * trn(ji,jj,jk,jpdoc)
223               ! tranfer of DOC to POC due to brownian motion
224               zaggdoc3 =  ( 5095. * trn(ji,jj,jk,jppoc) + 114. * 0.3 * trn(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trn(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     ! Total primary production per year
238     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,iksed+1) + sinking2(:,:,iksed+1) ) * e1e2t(:,:) * tmask(:,:,1) )
239     !
240     IF( ln_diatrc ) THEN
241         zrfact2 = 1.e3 * rfact2r
242         ik1  = iksed + 1
243         IF( lk_iomput ) THEN
244           IF( jnt == nrdttrc ) THEN
245              CALL iom_put( "EPC100"  , ( sinking(:,:,ik1) + sinking2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of carbon at 100m
246              CALL iom_put( "EPFE100" , ( sinkfer(:,:,ik1) + sinkfer2(:,:,ik1) ) * zrfact2 * tmask(:,:,1) ) ! Export of iron at 100m
247              CALL iom_put( "EPCAL100",   sinkcal(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of calcite  at 100m
248              CALL iom_put( "EPSI100" ,   sinksil(:,:,ik1)                       * zrfact2 * tmask(:,:,1) ) ! Export of biogenic silica at 100m
249           ENDIF
250         ELSE
251           trc2d(:,:,jp_pcs0_2d + 4) = sinking (:,:,ik1) * zrfact2 * tmask(:,:,1)
252           trc2d(:,:,jp_pcs0_2d + 5) = sinking2(:,:,ik1) * zrfact2 * tmask(:,:,1)
253           trc2d(:,:,jp_pcs0_2d + 6) = sinkfer (:,:,ik1) * zrfact2 * tmask(:,:,1)
254           trc2d(:,:,jp_pcs0_2d + 7) = sinkfer2(:,:,ik1) * zrfact2 * tmask(:,:,1)
255           trc2d(:,:,jp_pcs0_2d + 8) = sinksil (:,:,ik1) * zrfact2 * tmask(:,:,1)
256           trc2d(:,:,jp_pcs0_2d + 9) = sinkcal (:,:,ik1) * zrfact2 * tmask(:,:,1)
257         ENDIF
258      ENDIF
259      !
260      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
261         WRITE(charout, FMT="('sink')")
262         CALL prt_ctl_trc_info(charout)
263         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
264      ENDIF
265      !
266      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
267      !
268   END SUBROUTINE p4z_sink
269
270   SUBROUTINE p4z_sink_init
271      !!----------------------------------------------------------------------
272      !!                  ***  ROUTINE p4z_sink_init  ***
273      !!----------------------------------------------------------------------
274
275      t_oce_co2_exp = 0._wp
276      !
277   END SUBROUTINE p4z_sink_init
278
279#else
280   !!----------------------------------------------------------------------
281   !!   'Kriest sinking parameterisation'        key_kriest          ???
282   !!----------------------------------------------------------------------
283
284   SUBROUTINE p4z_sink ( kt, jnt )
285      !!---------------------------------------------------------------------
286      !!                ***  ROUTINE p4z_sink  ***
287      !!
288      !! ** Purpose :   Compute vertical flux of particulate matter due to
289      !!              gravitational sinking - Kriest parameterization
290      !!
291      !! ** Method  : - ???
292      !!---------------------------------------------------------------------
293      !
294      INTEGER, INTENT(in) :: kt, jnt
295      !
296      INTEGER  :: ji, jj, jk, jit, niter1, niter2
297      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh
298      REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc
299      REAL(wp) :: znum , zeps, zfm, zgm, zsm
300      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
301      REAL(wp) :: zval1, zval2, zval3, zval4
302      REAL(wp) :: zrfact2
303      INTEGER  :: ik1
304      CHARACTER (len=25) :: charout
305      REAL(wp), POINTER, DIMENSION(:,:,:) :: znum3d 
306      !!---------------------------------------------------------------------
307      !
308      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
309      !
310      CALL wrk_alloc( jpi, jpj, jpk, znum3d )
311      !
312      !     Initialisation of variables used to compute Sinking Speed
313      !     ---------------------------------------------------------
314
315      znum3d(:,:,:) = 0.e0
316      zval1 = 1. + xkr_zeta
317      zval2 = 1. + xkr_zeta + xkr_eta
318      zval3 = 1. + xkr_eta
319
320      !     Computation of the vertical sinking speed : Kriest et Evans, 2000
321      !     -----------------------------------------------------------------
322
323      DO jk = 1, jpkm1
324         DO jj = 1, jpj
325            DO ji = 1, jpi
326               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
327                  znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
328                  ! -------------- To avoid sinking speed over 50 m/day -------
329                  znum  = MIN( xnumm(jk), znum )
330                  znum  = MAX( 1.1      , znum )
331                  znum3d(ji,jj,jk) = znum
332                  !------------------------------------------------------------
333                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
334                  zfm   = xkr_frac**( 1. - zeps )
335                  zgm   = xkr_frac**( zval1 - zeps )
336                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
337                  zdiv1 = zeps - zval3
338                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
339                     &             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
340                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
341                     &             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
342                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
343               ENDIF
344            END DO
345         END DO
346      END DO
347
348      wscal(:,:,:) = MAX( wsbio3(:,:,:), 30._wp )
349
350      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
351      !   -----------------------------------------
352
353      sinking (:,:,:) = 0.e0
354      sinking2(:,:,:) = 0.e0
355      sinkcal (:,:,:) = 0.e0
356      sinkfer (:,:,:) = 0.e0
357      sinksil (:,:,:) = 0.e0
358
359     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
360     !   -----------------------------------------------------
361
362      niter1 = niter1max
363      niter2 = niter2max
364
365      DO jit = 1, niter1
366        CALL p4z_sink2( wsbio3, sinking , jppoc, niter1 )
367        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, niter1 )
368        CALL p4z_sink2( wscal , sinksil , jpgsi, niter1 )
369        CALL p4z_sink2( wscal , sinkcal , jpcal, niter1 )
370      END DO
371
372      DO jit = 1, niter2
373        CALL p4z_sink2( wsbio4, sinking2, jpnum, niter2 )
374      END DO
375
376     !  Exchange between organic matter compartments due to coagulation/disaggregation
377     !  ---------------------------------------------------
378
379      zval1 = 1. + xkr_zeta
380      zval2 = 1. + xkr_eta
381      zval3 = 3. + xkr_eta
382      zval4 = 4. + xkr_eta
383
384      DO jk = 1,jpkm1
385         DO jj = 1,jpj
386            DO ji = 1,jpi
387               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
388
389                  znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
390                  !-------------- To avoid sinking speed over 50 m/day -------
391                  znum  = min(xnumm(jk),znum)
392                  znum  = MAX( 1.1,znum)
393                  !------------------------------------------------------------
394                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
395                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
396                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
397                  zdiv2 = zeps - 2.
398                  zdiv3 = zeps - 3.
399                  zdiv4 = zeps - zval2
400                  zdiv5 = 2.* zeps - zval4
401                  zfm   = xkr_frac**( 1.- zeps )
402                  zsm   = xkr_frac**xkr_eta
403
404                  !    Part I : Coagulation dependant on turbulence
405                  !    ----------------------------------------------
406
407                  zagg1 =  0.163 * trn(ji,jj,jk,jpnum)**2               &
408                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
409                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
410                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
411                     &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
412                  zagg2 =  2*0.163*trn(ji,jj,jk,jpnum)**2*zfm*                       &
413                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
414                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
415                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
416                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
417                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
418                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))   
419
420                  zagg3 =  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3 
421                 
422                 !    Aggregation of small into large particles
423                 !    Part II : Differential settling
424                 !    ----------------------------------------------
425
426                  zagg4 =  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
427                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
428                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
429                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
430                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
431                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
432
433                  zagg5 =   2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2                         &
434                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
435                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
436                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
437                     &                 /zdiv) 
438
439                  !
440                  !     Fractionnation by swimming organisms
441                  !     ------------------------------------
442
443                  zfract = 2.*3.141*0.125*trn(ji,jj,jk,jpmes)*12./0.12/0.06**3*trn(ji,jj,jk,jpnum)  &
444                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
445                    &      * 10000.*xstep
446
447                  !     Aggregation of DOC to small particles
448                  !     --------------------------------------
449
450                  zaggdoc = 0.83 * trn(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)   &
451                     &        + 0.005 * 231. * trn(ji,jj,jk,jpdoc) * xstep * trn(ji,jj,jk,jpdoc)
452                  zaggdoc1 = 271. * trn(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  &
453                     &  + 0.02 * 16706. * trn(ji,jj,jk,jppoc) * xstep * trn(ji,jj,jk,jpdoc)
454
455# if defined key_degrad
456                   zagg1   = zagg1   * facvol(ji,jj,jk)                 
457                   zagg2   = zagg2   * facvol(ji,jj,jk)                 
458                   zagg3   = zagg3   * facvol(ji,jj,jk)                 
459                   zagg4   = zagg4   * facvol(ji,jj,jk)                 
460                   zagg5   = zagg5   * facvol(ji,jj,jk)                 
461                   zaggdoc = zaggdoc * facvol(ji,jj,jk)                 
462                   zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)
463# endif
464                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
465                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
466                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
467                  !
468                  znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
469                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
470                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
471                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1
472
473               ENDIF
474            END DO
475         END DO
476      END DO
477
478     ! Total primary production per year
479     t_oce_co2_exp = t_oce_co2_exp + glob_sum( ( sinking(:,:,:) ) * cvol(:,:,:) )
480     !
481      IF( ln_diatrc ) THEN
482         !
483         ik1 = iksed + 1
484         zrfact2 = 1.e3 * rfact2r
485         IF( jnt == nrdttrc ) THEN
486           CALL iom_put( "POCFlx"  , sinking (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! POC export
487           CALL iom_put( "NumFlx"  , sinking2 (:,:,:)     * zrfact2 * tmask(:,:,:) )  ! Num export
488           CALL iom_put( "SiFlx"   , sinksil (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Silica export
489           CALL iom_put( "CaCO3Flx", sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:) )  ! Calcite export
490           CALL iom_put( "xnum"    , znum3d  (:,:,:)                * tmask(:,:,:) )  ! Number of particles in aggregats
491           CALL iom_put( "W1"      , wsbio3  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of POC
492           CALL iom_put( "W2"      , wsbio4  (:,:,:)                * tmask(:,:,:) )  ! sinking speed of aggregats
493         ENDIF
494# if ! defined key_iomput
495         trc2d(:,:  ,jp_pcs0_2d + 4)  = sinking (:,:,ik1)    * zrfact2 * tmask(:,:,1)
496         trc2d(:,:  ,jp_pcs0_2d + 5)  = sinking2(:,:,ik1)    * zrfact2 * tmask(:,:,1)
497         trc2d(:,:  ,jp_pcs0_2d + 6)  = sinkfer (:,:,ik1)    * zrfact2 * tmask(:,:,1)
498         trc2d(:,:  ,jp_pcs0_2d + 7)  = sinksil (:,:,ik1)    * zrfact2 * tmask(:,:,1)
499         trc2d(:,:  ,jp_pcs0_2d + 8)  = sinkcal (:,:,ik1)    * zrfact2 * tmask(:,:,1)
500         trc3d(:,:,:,jp_pcs0_3d + 11) = sinking (:,:,:)      * zrfact2 * tmask(:,:,:)
501         trc3d(:,:,:,jp_pcs0_3d + 12) = sinking2(:,:,:)      * zrfact2 * tmask(:,:,:)
502         trc3d(:,:,:,jp_pcs0_3d + 13) = sinksil (:,:,:)      * zrfact2 * tmask(:,:,:)
503         trc3d(:,:,:,jp_pcs0_3d + 14) = sinkcal (:,:,:)      * zrfact2 * tmask(:,:,:)
504         trc3d(:,:,:,jp_pcs0_3d + 15) = znum3d  (:,:,:)                * tmask(:,:,:)
505         trc3d(:,:,:,jp_pcs0_3d + 16) = wsbio3  (:,:,:)                * tmask(:,:,:)
506         trc3d(:,:,:,jp_pcs0_3d + 17) = wsbio4  (:,:,:)                * tmask(:,:,:)
507# endif
508        !
509      ENDIF
510      !
511      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
512         WRITE(charout, FMT="('sink')")
513         CALL prt_ctl_trc_info(charout)
514         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
515      ENDIF
516      !
517      CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
518      !
519      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
520      !
521   END SUBROUTINE p4z_sink
522
523
524   SUBROUTINE p4z_sink_init
525      !!----------------------------------------------------------------------
526      !!                  ***  ROUTINE p4z_sink_init  ***
527      !!
528      !! ** Purpose :   Initialization of sinking parameters
529      !!                Kriest parameterization only
530      !!
531      !! ** Method  :   Read the nampiskrs namelist and check the parameters
532      !!      called at the first timestep
533      !!
534      !! ** input   :   Namelist nampiskrs
535      !!----------------------------------------------------------------------
536      INTEGER  ::   jk, jn, kiter
537      INTEGER  ::   ios                 ! Local integer output status for namelist read
538      REAL(wp) ::   znum, zdiv
539      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
540      REAL(wp) ::   zmin, zmax, zl, zr, xacc
541      !
542      NAMELIST/nampiskrs/ xkr_sfact, xkr_stick ,  &
543         &                xkr_nnano, xkr_ndiat, xkr_nmicro, xkr_nmeso, xkr_naggr
544      !!----------------------------------------------------------------------
545      !
546      IF( nn_timing == 1 )  CALL timing_start('p4z_sink_init')
547      !
548
549      REWIND( numnatp_ref )              ! Namelist nampiskrs in reference namelist : Pisces sinking Kriest
550      READ  ( numnatp_ref, nampiskrs, IOSTAT = ios, ERR = 901)
551901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in reference namelist', lwp )
552
553      REWIND( numnatp_cfg )              ! Namelist nampiskrs in configuration namelist : Pisces sinking Kriest
554      READ  ( numnatp_cfg, nampiskrs, IOSTAT = ios, ERR = 902 )
555902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampiskrs in configuration namelist', lwp )
556      IF(lwm) WRITE ( numonp, nampiskrs )
557
558      IF(lwp) THEN
559         WRITE(numout,*)
560         WRITE(numout,*) ' Namelist : nampiskrs'
561         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
562         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
563         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
564         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
565         WRITE(numout,*) '    Nbr of cell in microzoo size class       xkr_nmicro   = ', xkr_nmicro
566         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
567         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
568      ENDIF
569
570
571      ! max and min vertical particle speed
572      xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
573      xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
574      IF (lwp) WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
575
576      !
577      !    effect of the sizes of the different living pools on particle numbers
578      !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
579      !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
580      !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
581      !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
582      !    doc aggregates = 1um
583      ! ----------------------------------------------------------
584
585      xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
586      xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
587      xkr_dmicro = 1. / ( xkr_massp * xkr_nmicro )
588      xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
589      xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
590
591      !!---------------------------------------------------------------------
592      !!    'key_kriest'                                                  ???
593      !!---------------------------------------------------------------------
594      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
595      !  Search of the maximum number of particles in aggregates for each k-level.
596      !  Bissection Method
597      !--------------------------------------------------------------------
598      IF (lwp) THEN
599        WRITE(numout,*)
600        WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
601      ENDIF
602
603      xacc     =  0.001_wp
604      kiter    = 50
605      zmin     =  1.10_wp
606      zmax     = xkr_mass_max / xkr_mass_min
607      xkr_frac = zmax
608
609      DO jk = 1,jpk
610         zl = zmin
611         zr = zmax
612         wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2
613         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
614         znum = zl - 1.
615         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
616            & - ( xkr_wsbio_max * xkr_eta * znum * &
617            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
618            & - wmax
619
620         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
621         znum = zr - 1.
622         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
623            & - ( xkr_wsbio_max * xkr_eta * znum * &
624            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
625            & - wmax
626iflag:   DO jn = 1, kiter
627            IF    ( zwl == 0._wp ) THEN   ;   znummax = zl
628            ELSEIF( zwr == 0._wp ) THEN   ;   znummax = zr
629            ELSE
630               znummax = ( zr + zl ) / 2.
631               zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
632               znum = znummax - 1.
633               zws =  xkr_wsbio_min * xkr_zeta / zdiv &
634                  & - ( xkr_wsbio_max * xkr_eta * znum * &
635                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
636                  & - wmax
637               IF( zws * zwl < 0. ) THEN   ;   zr = znummax
638               ELSE                        ;   zl = znummax
639               ENDIF
640               zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
641               znum = zl - 1.
642               zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
643                  & - ( xkr_wsbio_max * xkr_eta * znum * &
644                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
645                  & - wmax
646
647               zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
648               znum = zr - 1.
649               zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
650                  & - ( xkr_wsbio_max * xkr_eta * znum * &
651                  &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
652                  & - wmax
653               !
654               IF ( ABS ( zws )  <= xacc ) EXIT iflag
655               !
656            ENDIF
657            !
658         END DO iflag
659
660         xnumm(jk) = znummax
661         IF (lwp) WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
662         !
663      END DO
664      !
665      t_oce_co2_exp = 0._wp
666      !
667      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink_init')
668      !
669  END SUBROUTINE p4z_sink_init
670
671#endif
672
673   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
674      !!---------------------------------------------------------------------
675      !!                     ***  ROUTINE p4z_sink2  ***
676      !!
677      !! ** Purpose :   Compute the sedimentation terms for the various sinking
678      !!     particles. The scheme used to compute the trends is based
679      !!     on MUSCL.
680      !!
681      !! ** Method  : - this ROUTINE compute not exactly the advection but the
682      !!      transport term, i.e.  div(u*tra).
683      !!---------------------------------------------------------------------
684      !
685      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
686      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting
687      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
688      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
689      !!
690      INTEGER  ::   ji, jj, jk, jn
691      REAL(wp) ::   zigma,zew,zign, zflx, zstep
692      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb 
693      !!---------------------------------------------------------------------
694      !
695      IF( nn_timing == 1 )  CALL timing_start('p4z_sink2')
696      !
697      ! Allocate temporary workspace
698      CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
699
700      zstep = rfact2 / FLOAT( kiter ) / 2.
701
702      ztraz(:,:,:) = 0.e0
703      zakz (:,:,:) = 0.e0
704      ztrb (:,:,:) = trn(:,:,:,jp_tra)
705
706      DO jk = 1, jpkm1
707         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
708      END DO
709      zwsink2(:,:,1) = 0.e0
710      IF( lk_degrad ) THEN
711         zwsink2(:,:,:) = zwsink2(:,:,:) * facvol(:,:,:)
712      ENDIF
713
714
715      ! Vertical advective flux
716      DO jn = 1, 2
717         !  first guess of the slopes interior values
718         DO jk = 2, jpkm1
719            ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
720         END DO
721         ztraz(:,:,1  ) = 0.0
722         ztraz(:,:,jpk) = 0.0
723
724         ! slopes
725         DO jk = 2, jpkm1
726            DO jj = 1,jpj
727               DO ji = 1, jpi
728                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
729                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
730               END DO
731            END DO
732         END DO
733         
734         ! Slopes limitation
735         DO jk = 2, jpkm1
736            DO jj = 1, jpj
737               DO ji = 1, jpi
738                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
739                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
740               END DO
741            END DO
742         END DO
743         
744         ! vertical advective flux
745         DO jk = 1, jpkm1
746            DO jj = 1, jpj     
747               DO ji = 1, jpi   
748                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
749                  zew   = zwsink2(ji,jj,jk+1)
750                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
751               END DO
752            END DO
753         END DO
754         !
755         ! Boundary conditions
756         psinkflx(:,:,1  ) = 0.e0
757         psinkflx(:,:,jpk) = 0.e0
758         
759         DO jk=1,jpkm1
760            DO jj = 1,jpj
761               DO ji = 1, jpi
762                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
763                  trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
764               END DO
765            END DO
766         END DO
767
768      ENDDO
769
770      DO jk = 1,jpkm1
771         DO jj = 1,jpj
772            DO ji = 1, jpi
773               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
774               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
775            END DO
776         END DO
777      END DO
778
779      trn(:,:,:,jp_tra) = ztrb(:,:,:)
780      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
781      !
782      CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
783      !
784      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink2')
785      !
786   END SUBROUTINE p4z_sink2
787
788
789   INTEGER FUNCTION p4z_sink_alloc()
790      !!----------------------------------------------------------------------
791      !!                     ***  ROUTINE p4z_sink_alloc  ***
792      !!----------------------------------------------------------------------
793      ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4  (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) ,     &
794         &      sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                      ,     &               
795         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                      ,     &               
796#if defined key_kriest
797         &      xnumm(jpk)                                                        ,     &               
798#else
799         &      sinkfer2(jpi,jpj,jpk)                                             ,     &               
800#endif
801         &      sinkfer(jpi,jpj,jpk)                                              , STAT=p4z_sink_alloc )               
802         !
803      IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
804      !
805   END FUNCTION p4z_sink_alloc
806   
807#else
808   !!======================================================================
809   !!  Dummy module :                                   No PISCES bio-model
810   !!======================================================================
811CONTAINS
812   SUBROUTINE p4z_sink                    ! Empty routine
813   END SUBROUTINE p4z_sink
814#endif 
815
816   !!======================================================================
817END MODULE  p4zsink
Note: See TracBrowser for help on using the repository browser.