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

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zagg.F90 @ 6841

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

Various bug fixes + explicit gamma function for lability
  • Property svn:executable set to *
File size: 14.7 KB
Line 
1MODULE p5zagg
2   !!======================================================================
3   !!                         ***  MODULE p5zagg  ***
4   !! TOP :  PISCES  aggregation of particles
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
9   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
10   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
11   !!----------------------------------------------------------------------
12#if defined key_pisces_quota
13   !!----------------------------------------------------------------------
14   !!   p5z_agg       :  Compute aggregation of particles
15   !!----------------------------------------------------------------------
16   USE oce_trc         !  shared variables between ocean and passive tracers
17   USE trc             !  passive tracers common variables
18   USE sms_pisces      !  PISCES Source Minus Sink variables
19   USE p5zsink         !  PISCES sinking of particles
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   p5z_agg         ! called in p5zbio.F90
28
29   !!* Substitution
30#  include "top_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
33   !! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38#if ! defined key_kriest
39   !!----------------------------------------------------------------------
40   !!   'standard particles parameterisation'                  ???
41   !!----------------------------------------------------------------------
42
43   SUBROUTINE p5z_agg ( kt, knt )
44      !!---------------------------------------------------------------------
45      !!                     ***  ROUTINE p5z_agg  ***
46      !!
47      !! ** Purpose :   Compute aggregation of particles
48      !!
49      !! ** Method  : - ???
50      !!---------------------------------------------------------------------
51      INTEGER, INTENT(in) :: kt, knt
52      INTEGER  ::   ji, jj, jk
53      REAL(wp) ::   zaggpoc1, zaggpoc2, zaggpoc3, zaggpoc4
54      REAL(wp) ::   zaggpoc , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
55      REAL(wp) ::   zaggpon , zaggdon, zaggdon2, zaggdon3
56      REAL(wp) ::   zaggpop, zaggdop, zaggdop2, zaggdop3
57      REAL(wp) ::   zaggtmp, zfact, zmax, zstep
58      CHARACTER (len=25) :: charout
59      !!---------------------------------------------------------------------
60      !
61      IF( nn_timing == 1 )  CALL timing_start('p5z_agg')
62      !
63      !  Exchange between organic matter compartments due to coagulation/disaggregation
64      !  ---------------------------------------------------
65      DO jk = 1, jpkm1
66         DO jj = 1, jpj
67            DO ji = 1, jpi
68               !
69               zstep = xstep 
70# if defined key_degrad
71               zstep = zstep * facvol(ji,jj,jk)
72# endif
73               zfact = zstep * xdiss(ji,jj,jk)
74               !  Part I : Coagulation dependent on turbulence
75               zaggtmp = 25.9  * zfact * trb(ji,jj,jk,jppoc)
76               zaggpoc1 = zaggtmp * trb(ji,jj,jk,jppoc)
77               zaggtmp = 4452. * zfact * trb(ji,jj,jk,jpgoc)
78               zaggpoc2 = zaggtmp * trb(ji,jj,jk,jppoc)
79
80               ! Part II : Differential settling
81
82               !  Aggregation of small into large particles
83               zaggtmp =  47.1 * zstep * trb(ji,jj,jk,jpgoc)
84               zaggpoc3 = zaggtmp * trb(ji,jj,jk,jppoc)
85               zaggtmp =  3.3  * zstep * trb(ji,jj,jk,jppoc)
86               zaggpoc4 = zaggtmp * trb(ji,jj,jk,jppoc)
87
88               zaggpoc   = zaggpoc1 + zaggpoc2 + zaggpoc3 + zaggpoc4
89               zaggpon = zaggpoc * trb(ji,jj,jk,jppon) / ( trb(ji,jj,jk,jppoc) + rtrn)
90               zaggpop = zaggpoc * trb(ji,jj,jk,jppop) / ( trb(ji,jj,jk,jppoc) + rtrn)
91               zaggfe = zaggpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc)  + rtrn )
92
93               ! Aggregation of DOC to POC :
94               ! 1st term is shear aggregation of DOC-DOC
95               ! 2nd term is shear aggregation of DOC-POC
96               ! 3rd term is differential settling of DOC-POC
97               zaggtmp = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
98               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) )
99               zaggdoc  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
100               zaggdon  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
101               zaggdop  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
102
103               ! transfer of DOC to GOC :
104               ! 1st term is shear aggregation
105               ! 2nd term is differential settling
106               zaggtmp = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc)
107               zaggdoc2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
108               zaggdon2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
109               zaggdop2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
110
111               ! tranfer of DOC to POC due to brownian motion
112               zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep
113               zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
114               zaggdon3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
115               zaggdop3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
116
117               !  Update the trends
118               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zaggpoc + zaggdoc + zaggdoc3
119               tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zaggpon + zaggdon + zaggdon3
120               tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zaggpop + zaggdop + zaggdop3
121               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zaggpoc + zaggdoc2
122               tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zaggpon + zaggdon2
123               tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + zaggpop + zaggdop2
124               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
125               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
126               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
127               tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zaggdon - zaggdon2 - zaggdon3
128               tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zaggdop - zaggdop2 - zaggdop3
129               !
130               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3
131               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zaggpoc + zaggdoc2
132               !
133            END DO
134         END DO
135      END DO
136      !
137      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
138         WRITE(charout, FMT="('agg')")
139         CALL prt_ctl_trc_info(charout)
140         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
141      ENDIF
142      !
143      IF( nn_timing == 1 )  CALL timing_stop('p5z_agg')
144      !
145   END SUBROUTINE p5z_agg
146
147#else
148   !!----------------------------------------------------------------------
149   !!   'Kriest parameterisation'        key_kriest          ???
150   !!----------------------------------------------------------------------
151
152   SUBROUTINE p5z_agg ( kt, knt )
153      !!---------------------------------------------------------------------
154      !!                ***  ROUTINE p5z_agg  ***
155      !!
156      !! ** Purpose :   Compute aggregation of particles
157      !!
158      !! ** Method  : - ???
159      !!---------------------------------------------------------------------
160      !
161      INTEGER, INTENT(in) :: kt, knt
162      !
163      INTEGER  :: ji, jj, jk
164      REAL(wp) :: zagg1, zagg2, zagg3, zagg4, zagg5, zfract, zaggsi, zaggsh
165      REAL(wp) :: zagg , zaggdoc, zaggdoc1, znumdoc
166      REAL(wp) :: znum , zeps, zfm, zgm, zsm, zfactn, zfactp
167      REAL(wp) :: zdiv , zdiv1, zdiv2, zdiv3, zdiv4, zdiv5
168      REAL(wp) :: zval1, zval2, zval3, zval4
169      CHARACTER (len=25) :: charout
170      !!---------------------------------------------------------------------
171      !
172      IF( nn_timing == 1 )  CALL timing_start('p5z_agg')
173      !
174      !  Exchange between organic matter compartments due to coagulation/disaggregation
175      !  ---------------------------------------------------
176
177      zval1 = 1. + xkr_zeta
178      zval2 = 1. + xkr_eta
179      zval3 = 3. + xkr_eta
180      zval4 = 4. + xkr_eta
181
182      DO jk = 1,jpkm1
183         DO jj = 1,jpj
184            DO ji = 1,jpi
185               IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
186
187                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp
188                  !-------------- To avoid sinking speed over 50 m/day -------
189                  znum  = min(xnumm(jk),znum)
190                  znum  = MAX( 1.1,znum)
191                  !------------------------------------------------------------
192                  zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
193                  zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
194                  zdiv1 = MAX( 1.e-4, ABS( zeps - 4.   ) ) * SIGN( 1., zeps - 4.    )
195                  zdiv2 = zeps - 2.
196                  zdiv3 = zeps - 3.
197                  zdiv4 = zeps - zval2
198                  zdiv5 = 2.* zeps - zval4
199                  zfm   = xkr_frac**( 1.- zeps )
200                  zsm   = xkr_frac**xkr_eta
201
202                  !    Part I : Coagulation dependant on turbulence
203                  !    ----------------------------------------------
204                  zagg1 =  0.163  * trb(ji,jj,jk,jpnum)**2                                 &
205                     &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
206                     &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
207                     &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
208                     &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
209                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm                        &
210                     &            * ((xkr_mass_max**3+3.*(xkr_mass_max**2            &
211                     &            * xkr_mass_min*(zeps-1.)/zdiv2                     &
212                     &            + xkr_mass_max*xkr_mass_min**2*(zeps-1.)/zdiv3)    &
213                     &            + xkr_mass_min**3*(zeps-1)/zdiv1)                  &
214                     &            - zfm*xkr_mass_max**3*(1.+3.*((zeps-1.)            &
215                     &            / (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))   
216
217                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3 
218                 
219                 !    Aggregation of small into large particles
220                 !    Part II : Differential settling
221                 !    ----------------------------------------------
222                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       &
223                     &                 xkr_wsbio_min*(zeps-1.)**2                         &
224                     &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
225                     &                 -(1.-zfm)/(zdiv*(zeps-1.)))-                       &
226                     &                 ((zfm*zfm*xkr_mass_max**2*zsm-xkr_mass_min**2)     &
227                     &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
228
229                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         &
230                     &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
231                     &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
232                     &                 /zdiv3-(xkr_mass_min**2-zfm*zsm*xkr_mass_max**2)    &
233                     &                 /zdiv) 
234
235                  !
236                  !     Fractionnation by swimming organisms
237                  !     ------------------------------------
238                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  &
239                    &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
240                    &      * 10000.*xstep
241
242                  !     Aggregation of DOC to small particles
243                  !     --------------------------------------
244                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   &
245                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)
246                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  &
247                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)
248
249# if defined key_degrad
250                   zagg1   = zagg1   * facvol(ji,jj,jk)                 
251                   zagg2   = zagg2   * facvol(ji,jj,jk)                 
252                   zagg3   = zagg3   * facvol(ji,jj,jk)                 
253                   zagg4   = zagg4   * facvol(ji,jj,jk)                 
254                   zagg5   = zagg5   * facvol(ji,jj,jk)                 
255                   zaggdoc = zaggdoc * facvol(ji,jj,jk)                 
256                   zaggdoc1 = zaggdoc1 * facvol(ji,jj,jk)
257# endif
258                  zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
259                  zaggsi = ( zagg4 + zagg5 ) * xstep / 10.
260                  zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
261                  !
262                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
263                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
264                  zfactn  = trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
265                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + ( zaggdoc + zaggdoc1 ) * zfactn
266                  zfactp  = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
267                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + ( zaggdoc + zaggdoc1 ) * zfactp
268                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
269                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc1
270
271               ENDIF
272            END DO
273         END DO
274      END DO
275     !
276     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
277         WRITE(charout, FMT="('agg')")
278         CALL prt_ctl_trc_info(charout)
279         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
280     ENDIF
281     !
282     CALL wrk_dealloc( jpi, jpj, jpk, znum3d )
283     !
284     IF( nn_timing == 1 )  CALL timing_stop('p5z_agg')
285     !
286   END SUBROUTINE p5z_agg
287
288#endif
289
290#else
291   !!======================================================================
292   !!  Dummy module :                                   No PISCES bio-model
293   !!======================================================================
294CONTAINS
295   SUBROUTINE p5z_agg                    ! Empty routine
296   END SUBROUTINE p5z_agg
297#endif 
298
299   !!======================================================================
300END MODULE  p5zagg
Note: See TracBrowser for help on using the repository browser.