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

source: trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90 @ 1329

Last change on this file since 1329 was 1329, checked in by cetlod, 15 years ago

update modules to take into account the mask land points in NetCDF outputs, see ticket:322
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 29.1 KB
Line 
1MODULE p4zsink
2   !!======================================================================
3   !!                         ***  MODULE p4zsink  ***
4   !! TOP :   PISCES Compute 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#if defined key_pisces
9   !!----------------------------------------------------------------------
10   !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
11   !!----------------------------------------------------------------------
12   USE trc
13   USE oce_trc         !
14   USE sms_pisces
15   USE prtctl_trc
16
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   p4z_sink    ! called in p4zbio.F90
22
23   !! * Shared module variables
24   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
25     wsbio3, wsbio4,      &    !: POC and GOC sinking speeds
26     wscal                     !: Calcite and BSi sinking speeds
27
28   !! * Module variables
29   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &   !:
30     sinking, sinking2,   &    !: POC sinking fluxes (different meanings depending on the parameterization
31     sinkcal, sinksil,    &    !: CaCO3 and BSi sinking fluxes
32     sinkfer                   !: Small BFe sinking flux
33
34   REAL(wp) ::   &
35     xstep , xstep2            !: Time step duration for biology
36
37#if  defined key_kriest
38   REAL(wp)          ::       &   
39      xkr_sfact    = 250.  ,  &   !: Sinking factor
40      xkr_stick    = 0.2   ,  &   !: Stickiness
41      xkr_nnano    = 2.337 ,  &   !: Nbr of cell in nano size class
42      xkr_ndiat    = 3.718 ,  &   !: Nbr of cell in diatoms size class
43      xkr_nmeso    = 7.147 ,  &   !: Nbr of cell in mesozoo  size class
44      xkr_naggr    = 9.877        !: Nbr of cell in aggregates  size class
45
46   REAL(wp)          ::       &   
47      xkr_frac
48
49   REAL(wp), PUBLIC ::        &
50      xkr_dnano            ,  &   !: Size of particles in nano pool
51      xkr_ddiat            ,  &   !: Size of particles in diatoms pool
52      xkr_dmeso            ,  &   !: Size of particles in mesozoo pool
53      xkr_daggr            ,  &   !: Size of particles in aggregates pool
54      xkr_wsbio_min        ,  &   !: min vertical particle speed
55      xkr_wsbio_max               !: max vertical particle speed
56
57   REAL(wp), PUBLIC, DIMENSION(jpk) ::   &   !:
58      xnumm                       !:     maximum number of particles in aggregates
59
60#endif
61
62#if ! defined key_kriest
63   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &   !:
64     sinkfer2                  !: Big Fe sinking flux
65#endif 
66
67   !!* Substitution
68#  include "domzgr_substitute.h90"
69   !!----------------------------------------------------------------------
70   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007)
71   !! $Id$
72   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74
75CONTAINS
76
77#if defined key_kriest
78
79   SUBROUTINE p4z_sink ( kt, jnt )
80      !!---------------------------------------------------------------------
81      !!                ***  ROUTINE p4z_sink  ***
82      !!
83      !! ** Purpose :   Compute vertical flux of particulate matter due to
84      !!              gravitational sinking - Kriest parameterization
85      !!
86      !! ** Method  : - ???
87      !!---------------------------------------------------------------------
88
89      INTEGER, INTENT(in) :: kt, jnt
90      INTEGER  :: ji, jj, jk
91      INTEGER  :: iksed
92      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zaggsi, zaggsh
93      REAL(wp) :: zagg , zaggdoc, znumdoc
94      REAL(wp) :: znum , zeps, zfm, zgm, zsm
95      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
96      REAL(wp) :: zval1, zval2, zval3, zval4
97#if defined key_trc_dia3d
98      REAL(wp) ::   zrfact2
99#endif
100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   znum3d
101      CHARACTER (len=25) :: charout
102
103      !!---------------------------------------------------------------------
104
105      IF( ( kt * jnt ) == nittrc000  )  THEN
106          CALL p4z_sink_init   ! Initialization (first time-step only)
107          xstep  = rfact2 / rjjss      ! Time step duration for biology
108          xstep2 = rfact2 /  2.
109      ENDIF
110
111!     Initialisation of variables used to compute Sinking Speed
112!     ---------------------------------------------------------
113
114       znum3d(:,:,:) = 0.e0
115       iksed = 10
116       zval1 = 1. + xkr_zeta
117       zval2 = 1. + xkr_zeta + xkr_eta
118       zval3 = 1. + xkr_eta
119
120!     Computation of the vertical sinking speed : Kriest et Evans, 2000
121!     -----------------------------------------------------------------
122
123      DO jk = 1, jpkm1
124         DO jj = 1, jpj
125            DO ji = 1, jpi
126               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
127                  znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
128! -------------- To avoid sinking speed over 50 m/day -------
129                  znum  = MIN( xnumm(jk), znum )
130                  znum  = MAX( 1.1      , znum )
131                  znum3d(ji,jj,jk) = znum
132!------------------------------------------------------------
133                  zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
134                  zfm   = xkr_frac**( 1. - zeps )
135                  zgm   = xkr_frac**( zval1 - zeps )
136                  zdiv  = MAX( 1.e-4, ABS( zeps - zval2 ) ) * SIGN( 1., ( zeps - zval2 ) )
137                  zdiv1 = zeps - zval3
138                  wsbio3(ji,jj,jk) = xkr_wsbio_min * ( zeps - zval1 ) / zdiv    &
139     &                             - xkr_wsbio_max *   zgm * xkr_eta  / zdiv
140                  wsbio4(ji,jj,jk) = xkr_wsbio_min *   ( zeps-1. )    / zdiv1   &
141     &                             - xkr_wsbio_max *   zfm * xkr_eta  / zdiv1
142                  IF( znum == 1.1)   wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
143               ENDIF
144            END DO
145         END DO
146      END DO
147
148      wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
149
150
151!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
152!   -----------------------------------------
153
154      sinking (:,:,:) = 0.e0
155      sinking2(:,:,:) = 0.e0
156      sinkcal (:,:,:) = 0.e0
157      sinkfer (:,:,:) = 0.e0
158      sinksil (:,:,:) = 0.e0
159
160!   Compute the sedimentation term using p4zsink2 for all
161!   the sinking particles
162!   -----------------------------------------------------
163
164      CALL p4z_sink2( wsbio3, sinking , jppoc )
165      CALL p4z_sink2( wsbio4, sinking2, jpnum )
166      CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
167      CALL p4z_sink2( wscal , sinksil , jpdsi )
168      CALL p4z_sink2( wscal , sinkcal , jpcal )
169
170!  Exchange between organic matter compartments due to
171!  coagulation/disaggregation
172!  ---------------------------------------------------
173
174      zval1 = 1. + xkr_zeta
175      zval2 = 1. + xkr_eta
176      zval3 = 3. + xkr_eta
177      zval4 = 4. + xkr_eta
178
179      DO jk = 1,jpkm1
180         DO jj = 1,jpj
181            DO ji = 1,jpi
182               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
183
184                  znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
185! -------------- To avoid sinking speed over 50 m/day -------
186                  znum  = min(xnumm(jk),znum)
187                  znum  = MAX( 1.1,znum)
188!------------------------------------------------------------
189                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
190                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
191                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
192                  zdiv2 = zeps - 2.
193                  zdiv3 = zeps - 3.
194                  zdiv4 = zeps - zval2
195                  zdiv5 = 2.* zeps - zval4
196                  zfm   = xkr_frac**( 1.- zeps )
197                  zsm   = xkr_frac**xkr_eta
198
199!    Part I : Coagulation dependant on turbulence
200!    ----------------------------------------------
201
202                  zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2               &
203                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
204                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
205                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
206                     &            * (zeps-1.)**2/(zdiv2*zdiv3))            &
207# if defined key_off_degrad
208                     &                 *facvol(ji,jj,jk)       &
209# endif
210                     &    )
211
212                  zagg2 = (  2*0.163*trn(ji,jj,jk,jpnum)**2*zfm*                       &
213                     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
214                     &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
215                     &                    +xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
216                     &                    +xkr_mass_min**3*(zeps-1)/zdiv1)                  &
217                     &                    -zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)/           &
218                     &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))      &
219#    if defined key_off_degrad
220                     &                 *facvol(ji,jj,jk)             &
221#    endif
222                     &    )
223
224                  zagg3 = (  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3   &
225#    if defined key_off_degrad
226                     &                 *facvol(ji,jj,jk)             &
227#    endif
228                     &    )
229
230                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
231
232!    Aggregation of small into large particles
233!    Part II : Differential settling
234!    ----------------------------------------------
235
236                  zagg4 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
237                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
238                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
239                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
240                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
241                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )                     &
242# if defined key_off_degrad
243                     &                 *facvol(ji,jj,jk)        &
244# endif
245                     &    )
246
247                  zagg5 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2                         &
248                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
249                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
250                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
251                     &                 /zdiv)                   &
252# if defined key_off_degrad
253                     &                 *facvol(ji,jj,jk)        &
254# endif
255                     &    )
256
257                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
258
259                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
260
261!     Aggregation of DOC to small particles
262!     --------------------------------------
263
264                  zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc)               &
265                     &        + 1018.  * trn(ji,jj,jk,jppoc)  ) * xstep    &
266# if defined key_off_degrad
267                     &        * facvol(ji,jj,jk)                              &
268# endif
269                     &        * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
270
271                  znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
272                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc
273                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zaggdoc * znumdoc - zagg
274                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc
275
276               ENDIF
277            END DO
278         END DO
279      END DO
280
281#if defined key_trc_diaadd
282      zrfact2 = 1.e3 * rfact2r
283      DO jj = 1, jpj 
284         DO ji = 1, jpi
285            trc2d(ji,jj, jp_pcs0_2d + 4) = sinking (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
286            trc2d(ji,jj, jp_pcs0_2d + 5) = sinking2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
287            trc2d(ji,jj, jp_pcs0_2d + 6) = sinkfer (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
288            trc2d(ji,jj, jp_pcs0_2d + 7) = sinksil (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
289            trc2d(ji,jj, jp_pcs0_2d + 8) = sinkcal (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
290         ENDDO
291      ENDDO
292#  if defined key_trc_dia3d
293      DO jk = 1, jpk
294         DO jj = 1, jpj
295            DO ji = 1, jpi
296               trc3d(ji,jj,jk,jp_pcs0_3d + 11) = sinking (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
297               trc3d(ji,jj,jk,jp_pcs0_3d + 12) = sinking2(ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
298               trc3d(ji,jj,jk,jp_pcs0_3d + 13) = sinksil (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
299               trc3d(ji,jj,jk,jp_pcs0_3d + 14) = sinkcal (ji,jj,jk) * zrfact2 * tmask(ji,jj,jk)
300               trc3d(ji,jj,jk,jp_pcs0_3d + 15) = znum3d  (ji,jj,jk)           * tmask(ji,jj,jk)
301               trc3d(ji,jj,jk,jp_pcs0_3d + 16) = wsbio3  (ji,jj,jk)           * tmask(ji,jj,jk)
302               trc3d(ji,jj,jk,jp_pcs0_3d + 17) = wsbio4  (ji,jj,jk)           * tmask(ji,jj,jk)
303            ENDDO
304         ENDDO
305      ENDDO
306#  endif
307
308#endif
309      !
310       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
311         WRITE(charout, FMT="('sink')")
312         CALL prt_ctl_trc_info(charout)
313         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
314       ENDIF
315
316   END SUBROUTINE p4z_sink
317
318   SUBROUTINE p4z_sink_init
319      !!----------------------------------------------------------------------
320      !!                  ***  ROUTINE p4z_sink_init  ***
321      !!
322      !! ** Purpose :   Initialization of sinking parameters
323      !!                Kriest parameterization only
324      !!
325      !! ** Method  :   Read the nampiskrs namelist and check the parameters
326      !!      called at the first timestep (nittrc000)
327      !!
328      !! ** input   :   Namelist nampiskrs
329      !!
330      !!----------------------------------------------------------------------
331      INTEGER  ::   jk, jn, kiter
332      REAL(wp) ::   znum, zdiv
333      REAL(wp) ::   zws, zwr, zwl,wmax, znummax
334      REAL(wp) ::   zmin, zmax, zl, zr, xacc
335
336      NAMELIST/nampiskrs/ xkr_sfact, xkr_stick ,  &
337         &                xkr_nnano, xkr_ndiat, xkr_nmeso, xkr_naggr
338
339      !!----------------------------------------------------------------------
340      REWIND( numnat )                     ! read nampiskrs
341      READ  ( numnat, nampiskrs )
342
343      IF(lwp) THEN
344         WRITE(numout,*)
345         WRITE(numout,*) ' Namelist : nampiskrs'
346         WRITE(numout,*) '    Sinking factor                           xkr_sfact    = ', xkr_sfact
347         WRITE(numout,*) '    Stickiness                               xkr_stick    = ', xkr_stick
348         WRITE(numout,*) '    Nbr of cell in nano size class           xkr_nnano    = ', xkr_nnano
349         WRITE(numout,*) '    Nbr of cell in diatoms size class        xkr_ndiat    = ', xkr_ndiat
350         WRITE(numout,*) '    Nbr of cell in mesozoo size class        xkr_nmeso    = ', xkr_nmeso
351         WRITE(numout,*) '    Nbr of cell in aggregates size class     xkr_naggr    = ', xkr_naggr
352     ENDIF
353
354
355     ! max and min vertical particle speed
356     xkr_wsbio_min = xkr_sfact * xkr_mass_min**xkr_eta
357     xkr_wsbio_max = xkr_sfact * xkr_mass_max**xkr_eta
358     WRITE(numout,*) ' max and min vertical particle speed ', xkr_wsbio_min, xkr_wsbio_max
359
360     !
361     !    effect of the sizes of the different living pools on particle numbers
362     !    nano = 2um-20um -> mean size=6.32 um -> ws=2.596 -> xnum=xnnano=2.337
363     !    diat and microzoo = 10um-200um -> 44.7 -> 8.732 -> xnum=xndiat=3.718
364     !    mesozoo = 200um-2mm -> 632.45 -> 45.14 -> xnum=xnmeso=7.147
365     !    aggregates = 200um-10mm -> 1414 -> 74.34 -> xnum=xnaggr=9.877
366     !    doc aggregates = 1um
367     ! ----------------------------------------------------------
368
369     xkr_dnano = 1. / ( xkr_massp * xkr_nnano )
370     xkr_ddiat = 1. / ( xkr_massp * xkr_ndiat )
371     xkr_dmeso = 1. / ( xkr_massp * xkr_nmeso )
372     xkr_daggr = 1. / ( xkr_massp * xkr_naggr )
373
374      !!---------------------------------------------------------------------
375      !!    'key_kriest'                                                  ???
376      !!---------------------------------------------------------------------
377      !  COMPUTATION OF THE VERTICAL PROFILE OF MAXIMUM SINKING SPEED
378      !  Search of the maximum number of particles in aggregates for each k-level.
379      !  Bissection Method
380      !--------------------------------------------------------------------
381      WRITE(numout,*)
382      WRITE(numout,*)'    kriest : Compute maximum number of particles in aggregates'
383
384      xacc     =  0.001
385      kiter    = 50
386      zmin     =  1.10
387      zmax     = xkr_mass_max / xkr_mass_min
388      xkr_frac = zmax
389
390      DO jk = 1,jpk
391         zl = zmin
392         zr = zmax
393         wmax = 0.5 * fse3t(1,1,jk) * rjjss / rfact2
394         zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
395         znum = zl - 1.
396         zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
397            & - ( xkr_wsbio_max * xkr_eta * znum * &
398            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
399            & - wmax
400
401         zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
402         znum = zr - 1.
403         zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
404            & - ( xkr_wsbio_max * xkr_eta * znum * &
405            &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
406            & - wmax
407iflag:  DO jn = 1, kiter
408           IF( zwl == 0.e0 ) THEN
409              znummax = zl
410           ELSE IF ( zwr == 0.e0 ) THEN
411              znummax = zr
412           ELSE
413              znummax = ( zr + zl ) / 2.
414              zdiv = xkr_zeta + xkr_eta - xkr_eta * znummax
415              znum = znummax - 1.
416              zws =  xkr_wsbio_min * xkr_zeta / zdiv &
417                 & - ( xkr_wsbio_max * xkr_eta * znum * &
418                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
419                 & - wmax
420              IF( zws * zwl < 0. ) THEN
421                 zr = znummax
422              ELSE
423                 zl = znummax
424              ENDIF
425              zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
426              znum = zl - 1.
427              zwl =  xkr_wsbio_min * xkr_zeta / zdiv &
428                 & - ( xkr_wsbio_max * xkr_eta * znum * &
429                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
430                 & - wmax
431
432              zdiv = xkr_zeta + xkr_eta - xkr_eta * zr
433              znum = zr - 1.
434              zwr =  xkr_wsbio_min * xkr_zeta / zdiv &
435                 & - ( xkr_wsbio_max * xkr_eta * znum * &
436                 &     xkr_frac**( -xkr_zeta / znum ) / zdiv ) &
437                 & - wmax
438
439              IF ( ABS ( zws )  <= xacc ) EXIT iflag
440
441           ENDIF
442
443        END DO iflag
444
445        xnumm(jk) = znummax
446        WRITE(numout,*) '       jk = ', jk, ' wmax = ', wmax,' xnum max = ', xnumm(jk)
447
448     END DO
449
450  END SUBROUTINE p4z_sink_init
451
452#else
453
454   SUBROUTINE p4z_sink ( kt, jnt )
455      !!---------------------------------------------------------------------
456      !!                     ***  ROUTINE p4z_sink  ***
457      !!
458      !! ** Purpose :   Compute vertical flux of particulate matter due to
459      !!                gravitational sinking
460      !!
461      !! ** Method  : - ???
462      !!---------------------------------------------------------------------
463      INTEGER, INTENT(in) :: kt, jnt
464      INTEGER  ::   ji, jj, jk
465      INTEGER  ::   iksed
466      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
467      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
468      REAL(wp) ::   zfact, zwsmax
469#if defined key_trc_dia3d
470      REAL(wp) ::   zrfact2
471#endif
472      CHARACTER (len=25) :: charout
473      !!---------------------------------------------------------------------
474
475      IF( ( kt * jnt ) == nittrc000  )  THEN
476          xstep  = rfact2 / rjjss      ! Timestep duration for biology
477          xstep2 = rfact2 /  2.
478      ENDIF
479
480!    Sinking speeds of detritus is increased with depth as shown
481!    by data and from the coagulation theory
482!    -----------------------------------------------------------
483
484      iksed = 10
485
486      DO jk = 1, jpkm1
487         DO jj = 1, jpj
488            DO ji=1,jpi
489               zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
490               wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
491            END DO
492         END DO
493      END DO
494
495!      LIMIT THE VALUES OF THE SINKING SPEEDS
496!      TO AVOID NUMERICAL INSTABILITIES
497
498      wsbio3(:,:,:) = wsbio
499!
500! OA Below, this is garbage. the ideal would be to find a time-splitting
501! OA algorithm that does not increase the computing cost by too much
502! OA In ROMS, I have included a time-splitting procedure. But it is
503! OA too expensive as the loop is computed globally. Thus, a small e3t
504! OA at one place determines the number of subtimesteps globally
505! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
506
507      DO jk = 1,jpkm1
508         DO jj = 1, jpj
509            DO ji = 1, jpi
510               zwsmax = 0.8 * fse3t(ji,jj,jk) / xstep
511               wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax )
512               wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax )
513            END DO
514         END DO
515      END DO
516
517      wscal(:,:,:) = wsbio4(:,:,:)
518
519!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
520!   -----------------------------------------
521
522      sinking (:,:,:) = 0.e0
523      sinking2(:,:,:) = 0.e0
524      sinkcal (:,:,:) = 0.e0
525      sinkfer (:,:,:) = 0.e0
526      sinksil (:,:,:) = 0.e0
527      sinkfer2(:,:,:) = 0.e0
528
529!   Compute the sedimentation term using p4zsink2 for all
530!   the sinking particles
531!   -----------------------------------------------------
532
533      CALL p4z_sink2( wsbio3, sinking , jppoc )
534      CALL p4z_sink2( wsbio3, sinkfer , jpsfe )
535      CALL p4z_sink2( wsbio4, sinking2, jpgoc )
536      CALL p4z_sink2( wsbio4, sinkfer2, jpbfe )
537      CALL p4z_sink2( wsbio4, sinksil , jpdsi )
538      CALL p4z_sink2( wscal , sinkcal , jpcal )
539
540!  Exchange between organic matter compartments due to
541!  coagulation/disaggregation
542!  ---------------------------------------------------
543
544      DO jk = 1, jpkm1
545         DO jj = 1, jpj
546            DO ji = 1, jpi
547               zfact = xstep * xdiss(ji,jj,jk)
548
549!    Part I : Coagulation dependent on turbulence
550!    ----------------------------------------------
551
552# if defined key_off_degrad
553               zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
554# else
555               zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
556# endif
557
558# if defined key_off_degrad
559               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
560# else
561               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
562# endif
563
564!    Aggregation of small into large particles
565!    Part II : Differential settling
566!    ----------------------------------------------
567
568# if defined key_off_degrad
569               zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
570# else
571               zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
572# endif
573
574# if defined key_off_degrad
575               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
576# else
577               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
578# endif
579
580               zagg   = zagg1 + zagg2 + zagg3 + zagg4
581               zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
582
583!     Aggregation of DOC to small particles
584!     --------------------------------------
585
586               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
587# if defined key_off_degrad
588                  &      * facvol(ji,jj,jk)                           &
589# endif
590                  &      * zfact * trn(ji,jj,jk,jpdoc)
591
592               zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
593# if defined key_off_degrad
594                  &        * facvol(ji,jj,jk)                            &
595# endif     
596                  &        * trn(ji,jj,jk,jpdoc)
597!
598!  Update the trends
599!  -----------------
600!
601               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
602               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
603               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
604               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
605               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2
606
607            END DO
608         END DO
609      END DO
610
611# if defined key_trc_diaadd
612      zrfact2 = 1.e3 * rfact2r
613      DO jj = 1, jpj
614         DO ji = 1, jpi
615            trc2d(ji,jj, jp_pcs0_2d + 4) = sinking (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
616            trc2d(ji,jj, jp_pcs0_2d + 5) = sinking2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
617            trc2d(ji,jj, jp_pcs0_2d + 6) = sinkfer (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
618            trc2d(ji,jj, jp_pcs0_2d + 7) = sinkfer2(ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
619            trc2d(ji,jj, jp_pcs0_2d + 8) = sinksil (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
620            trc2d(ji,jj, jp_pcs0_2d + 9) = sinkcal (ji,jj,iksed+1) * zrfact2 * tmask(ji,jj,1)
621         ENDDO
622      ENDDO
623# endif
624      !
625       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
626         WRITE(charout, FMT="('sink')")
627         CALL prt_ctl_trc_info(charout)
628         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
629       ENDIF
630
631   END SUBROUTINE p4z_sink
632
633#endif
634
635   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra )
636      !!---------------------------------------------------------------------
637      !!                     ***  ROUTINE p4z_sink2  ***
638      !!
639      !! ** Purpose :   Compute the sedimentation terms for the various sinking
640      !!     particles. The scheme used to compute the trends is based
641      !!     on MUSCL.
642      !!
643      !! ** Method  : - this ROUTINE compute not exactly the advection but the
644      !!      transport term, i.e.  div(u*tra).
645      !!---------------------------------------------------------------------
646      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
647      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
648      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
649      !!
650      INTEGER  ::   ji, jj, jk, jn
651      REAL(wp) ::   zigma,zew,zign, zflx
652      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztraz, zakz
653      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwsink2
654      !!---------------------------------------------------------------------
655
656
657      ztraz(:,:,:) = 0.e0
658      zakz (:,:,:) = 0.e0
659
660      DO jk = 1, jpkm1
661# if defined key_off_degrad
662         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1) * facvol(:,:,jk)
663# else
664         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1)
665# endif
666      END DO
667      zwsink2(:,:,1) = 0.e0
668
669
670      ! Vertical advective flux
671      DO jn = 1, 2
672         !  first guess of the slopes interior values
673         DO jk = 2, jpkm1
674            ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
675         END DO
676         ztraz(:,:,1  ) = 0.0
677         ztraz(:,:,jpk) = 0.0
678
679         ! slopes
680         DO jk = 2, jpkm1
681            DO jj = 1,jpj
682               DO ji = 1, jpi
683                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
684                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
685               END DO
686            END DO
687         END DO
688         
689         ! Slopes limitation
690         DO jk = 2, jpkm1
691            DO jj = 1, jpj
692               DO ji = 1, jpi
693                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
694                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
695               END DO
696            END DO
697         END DO
698         
699         ! vertical advective flux
700         DO jk = 1, jpkm1
701            DO jj = 1, jpj     
702               DO ji = 1, jpi   
703                  zigma = zwsink2(ji,jj,jk+1) * xstep2 / fse3w(ji,jj,jk+1)
704                  zew   = zwsink2(ji,jj,jk+1)
705                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * xstep2
706               END DO
707            END DO
708         END DO
709         !
710         ! Boundary conditions
711         psinkflx(:,:,1  ) = 0.e0
712         psinkflx(:,:,jpk) = 0.e0
713         
714         DO jk=1,jpkm1
715            DO jj = 1,jpj
716               DO ji = 1, jpi
717                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
718                  trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
719               END DO
720            END DO
721         END DO
722
723      ENDDO
724
725      DO jk=1,jpkm1
726         DO jj = 1,jpj
727            DO ji = 1, jpi
728               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
729               trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + 2. * zflx
730            END DO
731         END DO
732      END DO
733
734      trn(:,:,:,jp_tra) = trb(:,:,:,jp_tra)
735      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
736
737      !
738   END SUBROUTINE p4z_sink2
739
740#else
741   !!======================================================================
742   !!  Dummy module :                                   No PISCES bio-model
743   !!======================================================================
744CONTAINS
745   SUBROUTINE p4z_sink                    ! Empty routine
746   END SUBROUTINE p4z_sink
747#endif 
748
749   !!======================================================================
750END MODULE  p4zsink
Note: See TracBrowser for help on using the repository browser.