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 branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z – NEMO

source: branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90 @ 7068

Last change on this file since 7068 was 7068, checked in by cetlod, 7 years ago

ROBUST5_CNRS : implementation of part I of new TOP interface - 1st step -, see ticket #1782

File size: 17.5 KB
Line 
1MODULE p4zsink
2   !!======================================================================
3   !!                         ***  MODULE p4zsink  ***
4   !! TOP :  PISCES  vertical flux of particulate matter due to gravitational sinking
5   !!======================================================================
6   !! History :   1.0  !  2004     (O. Aumont) Original code
7   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
8   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
9   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
10   !!----------------------------------------------------------------------
11   !!   p4z_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
12   !!   p4z_sink_init  :  Unitialisation of sinking speed parameters
13   !!   p4z_sink_alloc :  Allocate sinking speed variables
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 prtctl_trc      !  print control for debugging
19   USE iom             !  I/O manager
20   USE lib_mpp
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   p4z_sink         ! called in p4zbio.F90
26   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
27   PUBLIC   p4z_sink_alloc
28
29   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio3   !: POC sinking speed
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wsbio4   !: GOC sinking speed
31   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   wscal    !: Calcite and BSi sinking speeds
32
33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinking, sinking2  !: POC sinking fluxes
34   !                                                          !  (different meanings depending on the parameterization)
35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkcal, sinksil   !: CaCO3 and BSi sinking fluxes
36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer            !: Small BFe sinking fluxes
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sinkfer2           !: Big iron sinking fluxes
38
39   INTEGER  :: ik100
40
41   !!----------------------------------------------------------------------
42   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
43   !! $Id: p4zsink.F90 3160 2011-11-20 14:27:18Z cetlod $
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   !!----------------------------------------------------------------------
49   !!   'standard sinking parameterisation'                  ???
50   !!----------------------------------------------------------------------
51
52   SUBROUTINE p4z_sink ( kt, knt )
53      !!---------------------------------------------------------------------
54      !!                     ***  ROUTINE p4z_sink  ***
55      !!
56      !! ** Purpose :   Compute vertical flux of particulate matter due to
57      !!                gravitational sinking
58      !!
59      !! ** Method  : - ???
60      !!---------------------------------------------------------------------
61      INTEGER, INTENT(in) :: kt, knt
62      INTEGER  ::   ji, jj, jk, jit
63      INTEGER  ::   iiter1, iiter2
64      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
65      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
66      REAL(wp) ::   zfact, zwsmax, zmax
67      CHARACTER (len=25) :: charout
68      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
69      REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
70      !!---------------------------------------------------------------------
71      !
72      IF( nn_timing == 1 )  CALL timing_start('p4z_sink')
73      !
74      !    Sinking speeds of detritus is increased with depth as shown
75      !    by data and from the coagulation theory
76      !    -----------------------------------------------------------
77      DO jk = 1, jpkm1
78         DO jj = 1, jpj
79            DO ji = 1,jpi
80               zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
81               zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / 5000._wp
82               wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
83            END DO
84         END DO
85      END DO
86
87      ! limit the values of the sinking speeds to avoid numerical instabilities 
88      wsbio3(:,:,:) = wsbio
89      wscal (:,:,:) = wsbio4(:,:,:)
90      !
91      ! OA This is (I hope) a temporary solution for the problem that may
92      ! OA arise in specific situation where the CFL criterion is broken
93      ! OA for vertical sedimentation of particles. To avoid this, a time
94      ! OA splitting algorithm has been coded. A specific maximum
95      ! OA iteration number is provided and may be specified in the namelist
96      ! OA This is to avoid very large iteration number when explicit free
97      ! OA surface is used (for instance). When niter?max is set to 1,
98      ! OA this computation is skipped. The crude old threshold method is
99      ! OA then applied. This also happens when niter exceeds nitermax.
100      IF( MAX( niter1max, niter2max ) == 1 ) THEN
101        iiter1 = 1
102        iiter2 = 1
103      ELSE
104        iiter1 = 1
105        iiter2 = 1
106        DO jk = 1, jpkm1
107          DO jj = 1, jpj
108             DO ji = 1, jpi
109                IF( tmask(ji,jj,jk) == 1) THEN
110                   zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep
111                   iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
112                   iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
113                ENDIF
114             END DO
115          END DO
116        END DO
117        IF( lk_mpp ) THEN
118           CALL mpp_max( iiter1 )
119           CALL mpp_max( iiter2 )
120        ENDIF
121        iiter1 = MIN( iiter1, niter1max )
122        iiter2 = MIN( iiter2, niter2max )
123      ENDIF
124
125      DO jk = 1,jpkm1
126         DO jj = 1, jpj
127            DO ji = 1, jpi
128               IF( tmask(ji,jj,jk) == 1 ) THEN
129                 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
130                 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
131                 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
132               ENDIF
133            END DO
134         END DO
135      END DO
136
137      !  Initializa to zero all the sinking arrays
138      !   -----------------------------------------
139      sinking (:,:,:) = 0.e0
140      sinking2(:,:,:) = 0.e0
141      sinkcal (:,:,:) = 0.e0
142      sinkfer (:,:,:) = 0.e0
143      sinksil (:,:,:) = 0.e0
144      sinkfer2(:,:,:) = 0.e0
145
146      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
147      !   -----------------------------------------------------
148      DO jit = 1, iiter1
149        CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
150        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
151      END DO
152
153      DO jit = 1, iiter2
154        CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
155        CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
156        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
157        CALL p4z_sink2( wscal , sinkcal , jpcal, iiter2 )
158      END DO
159
160      !  Exchange between organic matter compartments due to coagulation/disaggregation
161      !  ---------------------------------------------------
162      DO jk = 1, jpkm1
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               !
166               zfact = xstep * xdiss(ji,jj,jk)
167               !  Part I : Coagulation dependent on turbulence
168               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
169               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
170
171               ! Part II : Differential settling
172
173               !  Aggregation of small into large particles
174               zagg3 =  47.1 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
175               zagg4 =  3.3  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
176
177               zagg   = zagg1 + zagg2 + zagg3 + zagg4
178               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
179
180               ! Aggregation of DOC to POC :
181               ! 1st term is shear aggregation of DOC-DOC
182               ! 2nd term is shear aggregation of DOC-POC
183               ! 3rd term is differential settling of DOC-POC
184               zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
185               &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
186               ! transfer of DOC to GOC :
187               ! 1st term is shear aggregation
188               ! 2nd term is differential settling
189               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
190               ! tranfer of DOC to POC due to brownian motion
191               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) * xstep * 0.3 * trb(ji,jj,jk,jpdoc)
192
193               !  Update the trends
194               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
195               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
196               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
197               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
198               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
199               !
200            END DO
201         END DO
202      END DO
203
204
205     ! Total carbon export per year
206     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
207        &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
208     !
209     IF( lk_iomput ) THEN
210       IF( knt == nrdttrc ) THEN
211          CALL wrk_alloc( jpi, jpj,      zw2d )
212          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
213          zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
214          !
215          IF( iom_use( "EPC100" ) )  THEN
216              zw2d(:,:) = ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of carbon at 100m
217              CALL iom_put( "EPC100"  , zw2d )
218          ENDIF
219          IF( iom_use( "EPFE100" ) )  THEN
220              zw2d(:,:) = ( sinkfer(:,:,ik100) + sinkfer2(:,:,ik100) ) * zfact * tmask(:,:,1) ! Export of iron at 100m
221              CALL iom_put( "EPFE100"  , zw2d )
222          ENDIF
223          IF( iom_use( "EPCAL100" ) )  THEN
224              zw2d(:,:) = sinkcal(:,:,ik100) * zfact * tmask(:,:,1) ! Export of calcite at 100m
225              CALL iom_put( "EPCAL100"  , zw2d )
226          ENDIF
227          IF( iom_use( "EPSI100" ) )  THEN
228              zw2d(:,:) =  sinksil(:,:,ik100) * zfact * tmask(:,:,1) ! Export of bigenic silica at 100m
229              CALL iom_put( "EPSI100"  , zw2d )
230          ENDIF
231          IF( iom_use( "EXPC" ) )  THEN
232              zw3d(:,:,:) = ( sinking(:,:,:) + sinking2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of carbon in the water column
233              CALL iom_put( "EXPC"  , zw3d )
234          ENDIF
235          IF( iom_use( "EXPFE" ) )  THEN
236              zw3d(:,:,:) = ( sinkfer(:,:,:) + sinkfer2(:,:,:) ) * zfact * tmask(:,:,:) ! Export of iron
237              CALL iom_put( "EXPFE"  , zw3d )
238          ENDIF
239          IF( iom_use( "EXPCAL" ) )  THEN
240              zw3d(:,:,:) = sinkcal(:,:,:) * zfact * tmask(:,:,:) ! Export of calcite
241              CALL iom_put( "EXPCAL"  , zw3d )
242          ENDIF
243          IF( iom_use( "EXPSI" ) )  THEN
244              zw3d(:,:,:) = sinksil(:,:,:) * zfact * tmask(:,:,:) ! Export of bigenic silica
245              CALL iom_put( "EXPSI"  , zw3d )
246          ENDIF
247          IF( iom_use( "tcexp" ) )  CALL iom_put( "tcexp" , t_oce_co2_exp * zfact )   ! molC/s
248          !
249          CALL wrk_dealloc( jpi, jpj,      zw2d )
250          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
251        ENDIF
252      ENDIF
253      !
254      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
255         WRITE(charout, FMT="('sink')")
256         CALL prt_ctl_trc_info(charout)
257         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
258      ENDIF
259      !
260      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink')
261      !
262   END SUBROUTINE p4z_sink
263
264   SUBROUTINE p4z_sink_init
265      !!----------------------------------------------------------------------
266      !!                  ***  ROUTINE p4z_sink_init  ***
267      !!----------------------------------------------------------------------
268      INTEGER :: jk
269
270      ik100 = 10        !  last level where depth less than 100 m
271      DO jk = jpkm1, 1, -1
272         IF( gdept_1d(jk) > 100. )  ik100 = jk - 1
273      END DO
274      IF (lwp) WRITE(numout,*)
275      IF (lwp) WRITE(numout,*) ' Level corresponding to 100m depth ',  ik100 + 1
276      IF (lwp) WRITE(numout,*)
277      !
278      t_oce_co2_exp = 0._wp
279      !
280   END SUBROUTINE p4z_sink_init
281
282   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
283      !!---------------------------------------------------------------------
284      !!                     ***  ROUTINE p4z_sink2  ***
285      !!
286      !! ** Purpose :   Compute the sedimentation terms for the various sinking
287      !!     particles. The scheme used to compute the trends is based
288      !!     on MUSCL.
289      !!
290      !! ** Method  : - this ROUTINE compute not exactly the advection but the
291      !!      transport term, i.e.  div(u*tra).
292      !!---------------------------------------------------------------------
293      !
294      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index     
295      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting
296      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
297      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
298      !!
299      INTEGER  ::   ji, jj, jk, jn
300      REAL(wp) ::   zigma,zew,zign, zflx, zstep
301      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztraz, zakz, zwsink2, ztrb 
302      !!---------------------------------------------------------------------
303      !
304      IF( nn_timing == 1 )  CALL timing_start('p4z_sink2')
305      !
306      ! Allocate temporary workspace
307      CALL wrk_alloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
308
309      zstep = rfact2 / FLOAT( kiter ) / 2.
310
311      ztraz(:,:,:) = 0.e0
312      zakz (:,:,:) = 0.e0
313      ztrb (:,:,:) = trb(:,:,:,jp_tra)
314
315      DO jk = 1, jpkm1
316         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
317      END DO
318      zwsink2(:,:,1) = 0.e0
319
320
321      ! Vertical advective flux
322      DO jn = 1, 2
323         !  first guess of the slopes interior values
324         DO jk = 2, jpkm1
325            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
326         END DO
327         ztraz(:,:,1  ) = 0.0
328         ztraz(:,:,jpk) = 0.0
329
330         ! slopes
331         DO jk = 2, jpkm1
332            DO jj = 1,jpj
333               DO ji = 1, jpi
334                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
335                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
336               END DO
337            END DO
338         END DO
339         
340         ! Slopes limitation
341         DO jk = 2, jpkm1
342            DO jj = 1, jpj
343               DO ji = 1, jpi
344                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
345                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
346               END DO
347            END DO
348         END DO
349         
350         ! vertical advective flux
351         DO jk = 1, jpkm1
352            DO jj = 1, jpj     
353               DO ji = 1, jpi   
354                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
355                  zew   = zwsink2(ji,jj,jk+1)
356                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
357               END DO
358            END DO
359         END DO
360         !
361         ! Boundary conditions
362         psinkflx(:,:,1  ) = 0.e0
363         psinkflx(:,:,jpk) = 0.e0
364         
365         DO jk=1,jpkm1
366            DO jj = 1,jpj
367               DO ji = 1, jpi
368                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
369                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
370               END DO
371            END DO
372         END DO
373
374      ENDDO
375
376      DO jk = 1,jpkm1
377         DO jj = 1,jpj
378            DO ji = 1, jpi
379               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
380               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
381            END DO
382         END DO
383      END DO
384
385      trb(:,:,:,jp_tra) = ztrb(:,:,:)
386      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
387      !
388      CALL wrk_dealloc( jpi, jpj, jpk, ztraz, zakz, zwsink2, ztrb )
389      !
390      IF( nn_timing == 1 )  CALL timing_stop('p4z_sink2')
391      !
392   END SUBROUTINE p4z_sink2
393
394
395   INTEGER FUNCTION p4z_sink_alloc()
396      !!----------------------------------------------------------------------
397      !!                     ***  ROUTINE p4z_sink_alloc  ***
398      !!----------------------------------------------------------------------
399      ALLOCATE( wsbio3 (jpi,jpj,jpk) , wsbio4  (jpi,jpj,jpk) , wscal(jpi,jpj,jpk) ,     &
400         &      sinking(jpi,jpj,jpk) , sinking2(jpi,jpj,jpk)                      ,     &               
401         &      sinkcal(jpi,jpj,jpk) , sinksil (jpi,jpj,jpk)                      ,     &               
402         &      sinkfer2(jpi,jpj,jpk)                                             ,     &               
403         &      sinkfer(jpi,jpj,jpk)                                              , STAT=p4z_sink_alloc )               
404         !
405      IF( p4z_sink_alloc /= 0 ) CALL ctl_warn('p4z_sink_alloc : failed to allocate arrays.')
406      !
407   END FUNCTION p4z_sink_alloc
408   
409   !!======================================================================
410END MODULE p4zsink
Note: See TracBrowser for help on using the repository browser.