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_r12563_ASINTER-06_ABL_improvement/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/TOP/TRP/trcsink.F90 @ 12808

Last change on this file since 12808 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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, 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., 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.