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

source: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90 @ 7180

Last change on this file since 7180 was 7180, checked in by aumont, 7 years ago

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