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.
trcsink.F90 in NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/TOP/TRP/trcsink.F90 @ 13247

Last change on this file since 13247 was 13247, checked in by francesca, 4 years ago

dev_r12558_HPC-08_epico_Extra_Halo: merge with trunk@13227, see #2366
File size: 9.8 KB
Line 
1MODULE trcsink
2   !!======================================================================
3   !!                         ***  MODULE trcsink  ***
4   !! TOP :  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   !!             4.0  !  2018-12  (O. Aumont) Generalize the PISCES code to make it usable by any model
11   !!----------------------------------------------------------------------
12   !!   trc_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
13   !!----------------------------------------------------------------------
14   USE oce_trc         !  shared variables between ocean and passive tracers
15   USE trc             !  passive tracers common variables
16   USE lib_mpp
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC trc_sink
22   PUBLIC trc_sink_ini
23
24   INTEGER, PUBLIC :: nitermax      !: Maximum number of iterations for sinking
25
26   !! * Substitutions
27#  include "do_loop_substitute.h90"
28   !!----------------------------------------------------------------------
29   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
30   !! $Id: trcsink.F90 10069 2018-08-28 14:12:24Z nicolasmartin $
31   !! Software governed by the CeCILL license (see ./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   !!----------------------------------------------------------------------
36   !!   'standard sinking parameterisation'                  ???
37   !!----------------------------------------------------------------------
38
39   SUBROUTINE trc_sink ( kt, Kbb, Kmm, pwsink, psinkflx, jp_tra, rsfact )
40      !!---------------------------------------------------------------------
41      !!                     ***  ROUTINE trc_sink  ***
42      !!
43      !! ** Purpose :   Compute vertical flux of particulate matter due to
44      !!                gravitational sinking
45      !!
46      !! ** Method  : - ???
47      !!---------------------------------------------------------------------
48      INTEGER , INTENT(in)  :: kt
49      INTEGER , INTENT(in)  :: Kbb, Kmm
50      INTEGER , INTENT(in)  :: jp_tra    ! tracer index index     
51      REAL(wp), INTENT(in)  :: rsfact    ! time step duration
52      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pwsink
53      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx
54      INTEGER  ::   ji, jj, jk
55      INTEGER, DIMENSION(jpi, jpj) ::   iiter
56      REAL(wp) ::   zfact, zwsmax, zmax
57      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink
58      !!---------------------------------------------------------------------
59      !
60      IF( ln_timing )   CALL timing_start('trc_sink')
61      !
62      !
63      ! OA This is (I hope) a temporary solution for the problem that may
64      ! OA arise in specific situation where the CFL criterion is broken
65      ! OA for vertical sedimentation of particles. To avoid this, a time
66      ! OA splitting algorithm has been coded. A specific maximum
67      ! OA iteration number is provided and may be specified in the namelist
68      ! OA This is to avoid very large iteration number when explicit free
69      ! OA surface is used (for instance). When niter?max is set to 1,
70      ! OA this computation is skipped. The crude old threshold method is
71      ! OA then applied. This also happens when niter exceeds nitermax.
72      IF( nitermax == 1 ) THEN
73         iiter(:,:) = 1
74      ELSE
75         DO_2D_11_11
76            iiter(ji,jj) = 1
77            DO jk = 1, jpkm1
78               IF( tmask(ji,jj,jk) == 1.0 ) THEN
79                   zwsmax =  0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
80                   iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) )
81               ENDIF
82            END DO
83         END_2D
84         iiter(:,:) = MIN( iiter(:,:), nitermax )
85      ENDIF
86
87      DO_3D_11_11( 1,jpkm1 )
88         IF( tmask(ji,jj,jk) == 1.0 ) THEN
89           zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
90           zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
91         ELSE
92           ! provide a default value so there is no use of undefinite value in trc_sink2 for zwsink2 initialization
93           zwsink(ji,jj,jk) = 0.
94         ENDIF
95      END_3D
96
97      !  Initializa to zero all the sinking arrays
98      !  -----------------------------------------
99      psinkflx(:,:,:) = 0.e0
100
101      !   Compute the sedimentation term using trc_sink2 for the considered sinking particle
102      !   -----------------------------------------------------
103      CALL trc_sink2( Kbb, Kmm, zwsink, psinkflx, jp_tra, iiter, rsfact )
104      !
105      IF( ln_timing )   CALL timing_stop('trc_sink')
106      !
107   END SUBROUTINE trc_sink
108
109   SUBROUTINE trc_sink2( Kbb, Kmm, pwsink, psinkflx, jp_tra, kiter, rsfact )
110      !!---------------------------------------------------------------------
111      !!                     ***  ROUTINE trc_sink2  ***
112      !!
113      !! ** Purpose :   Compute the sedimentation terms for the various sinking
114      !!     particles. The scheme used to compute the trends is based
115      !!     on MUSCL.
116      !!
117      !! ** Method  : - this ROUTINE compute not exactly the advection but the
118      !!      transport term, i.e.  div(u*tra).
119      !!---------------------------------------------------------------------
120      INTEGER,  INTENT(in   )                         ::   Kbb, Kmm  ! time level indices
121      INTEGER,  INTENT(in   )                         ::   jp_tra    ! tracer index index     
122      REAL(wp), INTENT(in   )                         ::   rsfact    ! duration of time step
123      INTEGER,  INTENT(in   ), DIMENSION(jpi,jpj)     ::   kiter     ! number of iterations for time-splitting
124      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
125      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
126      !
127      INTEGER  ::   ji, jj, jk, jn
128      REAL(wp) ::   zigma,zew,zign, zflx, zstep
129      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb 
130      !!---------------------------------------------------------------------
131      !
132      IF( ln_timing )   CALL timing_start('trc_sink2')
133      !
134      ztraz(:,:,:) = 0.e0
135      zakz (:,:,:) = 0.e0
136      ztrb (:,:,:) = tr(:,:,:,jp_tra,Kbb)
137
138      DO jk = 1, jpkm1
139         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
140      END DO
141      zwsink2(:,:,1) = 0.e0
142
143
144      ! Vertical advective flux
145      DO jn = 1, 2
146         !  first guess of the slopes interior values
147         DO_2D_11_11
148            !
149            zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
150            !             
151            DO jk = 2, jpkm1
152               ztraz(ji,jj,jk) = ( tr(ji,jj,jk-1,jp_tra,Kbb) - tr(ji,jj,jk,jp_tra,Kbb) ) * tmask(ji,jj,jk)
153            END DO
154            ztraz(ji,jj,1  ) = 0.0
155            ztraz(ji,jj,jpk) = 0.0
156
157            ! slopes
158            DO jk = 2, jpkm1
159               zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
160               zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
161            END DO
162     
163            ! Slopes limitation
164            DO jk = 2, jpkm1
165               zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        &
166                  &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
167            END DO
168     
169            ! vertical advective flux
170            DO jk = 1, jpkm1
171               zigma = zwsink2(ji,jj,jk+1) * zstep / e3w(ji,jj,jk+1,Kmm)
172               zew   = zwsink2(ji,jj,jk+1)
173               psinkflx(ji,jj,jk+1) = -zew * ( tr(ji,jj,jk,jp_tra,Kbb) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
174            END DO
175            !
176            ! Boundary conditions
177            psinkflx(ji,jj,1  ) = 0.e0
178            psinkflx(ji,jj,jpk) = 0.e0
179     
180            DO jk=1,jpkm1
181               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
182               tr(ji,jj,jk,jp_tra,Kbb) = tr(ji,jj,jk,jp_tra,Kbb) + zflx
183            END DO
184         END_2D
185      END DO
186
187      DO_3D_11_11( 1,jpkm1 )
188         zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
189         ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
190      END_3D
191
192      tr(:,:,:,jp_tra,Kbb) = ztrb(:,:,:)
193      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
194      !
195      IF( ln_timing )  CALL timing_stop('trc_sink2')
196      !
197   END SUBROUTINE trc_sink2
198
199  SUBROUTINE trc_sink_ini
200      !!---------------------------------------------------------------------
201      !!                  ***  ROUTINE trc_sink_ini ***
202      !!
203      !! ** Purpose :   read  namelist options
204      !!----------------------------------------------------------------------
205      INTEGER ::   ios   ! Local integer output status for namelist read
206      !!
207      NAMELIST/namtrc_snk/ nitermax
208      !!----------------------------------------------------------------------
209      !
210      READ  ( numnat_ref, namtrc_snk, IOSTAT = ios, ERR = 907)
211907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_snk in reference namelist' )
212      READ  ( numnat_cfg, namtrc_snk, IOSTAT = ios, ERR = 908 )
213908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_snk in configuration namelist' )
214      IF(lwm) WRITE( numont, namtrc_snk )
215
216      IF(lwp) THEN                     !   ! Control print
217         WRITE(numout,*)
218         WRITE(numout,*) 'trc_sink : Sedimentation of particles '
219         WRITE(numout,*) '~~~~~~~ '
220         WRITE(numout,*) '   Namelist namtrc_snk : sedimentation of particles'
221         WRITE(numout,*) '   Maximum number of iterations   nitermax = ', nitermax
222         WRITE(numout,*)
223      ENDIF
224
225   END SUBROUTINE trc_sink_ini
226
227   !!======================================================================
228END MODULE trcsink
Note: See TracBrowser for help on using the repository browser.