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_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/TOP/TRP/trcsink.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
File size: 10.0 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      READ  ( numnat_ref, namtrc_snk, IOSTAT = ios, ERR = 907)
219907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_snk in reference namelist' )
220      READ  ( numnat_cfg, namtrc_snk, IOSTAT = ios, ERR = 908 )
221908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_snk in configuration namelist' )
222      IF(lwm) WRITE( numont, namtrc_snk )
223
224      IF(lwp) THEN                     !   ! Control print
225         WRITE(numout,*)
226         WRITE(numout,*) 'trc_sink : Sedimentation of particles '
227         WRITE(numout,*) '~~~~~~~ '
228         WRITE(numout,*) '   Namelist namtrc_snk : sedimentation of particles'
229         WRITE(numout,*) '   Maximum number of iterations   nitermax = ', nitermax
230         WRITE(numout,*)
231      ENDIF
232
233   END SUBROUTINE trc_sink_ini
234
235   !!======================================================================
236END MODULE trcsink
Note: See TracBrowser for help on using the repository browser.