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

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

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