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/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/TOP/TRP/trcsink.F90 @ 11573

Last change on this file since 11573 was 11573, checked in by jchanut, 5 years ago

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