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

Last change on this file since 7041 was 7041, checked in by cetlod, 4 years ago

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

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