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.
trczdf_imp.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trczdf_imp.F90 @ 650

Last change on this file since 650 was 615, checked in by opalod, 17 years ago

nemo_v2_bugfix_023:CE:change key_trc_ldfiso to key_trcldf_iso
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.9 KB
Line 
1MODULE trczdf_imp
2   !!==============================================================================
3   !!                    ***  MODULE  trczdf_imp  ***
4   !! Ocean passive tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_zdf_imp  : update the tracer trend with the vertical diffusion
9   !!                  using an implicit time-stepping.
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce_trc             ! ocean dynamics and active tracers variables
13   USE trc                 ! ocean passive tracers variables
14   USE trctrp_lec          ! passive tracers transport
15   USE prtctl_trc
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Routine accessibility
21   PUBLIC trc_zdf_imp          ! routine called by step.F90
22
23   !! * Module variable
24   REAL(wp), DIMENSION(jpk) ::   &
25      rdttrc                     ! vertical profile of 2 x tracer time-step
26
27   !! * Substitutions
28#  include "passivetrc_substitute.h90"
29   !!----------------------------------------------------------------------
30   !!   TOP 1.0 , LOCEAN-IPSL (2005)
31   !! $Header$
32   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
33   !!----------------------------------------------------------------------
34
35CONTAINS
36
37   SUBROUTINE trc_zdf_imp( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE trc_zdf_imp  ***
40      !!
41      !! ** Purpose :   Compute the trend due to the vertical tracer mixing
42      !!      using an implicit time stepping and add it to the general trend
43      !!      of the tracer equations.
44      !!
45      !! ** Method  :   The vertical diffusion of tracers tra is given by:
46      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(tra) )
47      !!      It is thus evaluated using a backward time scheme
48      !!      Surface and bottom boundary conditions: no diffusive flux on
49      !!      both tracers (bottom, applied through the masked field avt).
50      !!      Add this trend to the general trend tra :
51      !!          tra = tra + dz( avt dz(t) )
52      !!         (tra = tra + dz( avs dz(t) ) if lk_zdfddmtrc=T)
53      !!
54      !! ** Action  : - Update tra with the before vertical diffusion trend
55      !!              - save the trends in trtrd ('key_trc_diatrd')
56      !!
57      !! History :
58      !!   6.0  !  90-10  (B. Blanke)  Original code
59      !!   7.0  !  91-11  (G. Madec)
60      !!        !  92-06  (M. Imbard)  correction on tracer trend loops
61      !!        !  96-01  (G. Madec)  statement function for e3
62      !!        !  97-05  (G. Madec)  vertical component of isopycnal
63      !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
64      !!        !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
65      !!        !  00-05  (MA Foujols) add lbc for tracer trends
66      !!        !  00-06  (O Aumont)  correct isopycnal scheme suppress
67      !!        !                     avt multiple correction
68      !!        !  00-08  (G. Madec)  double diffusive mixing
69      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
70      !!   9.0  !  04-03  (C. Ethe )  adapted for passive tracers
71      !!---------------------------------------------------------------------
72      !! * Arguments
73      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
74      INTEGER ::   ikst, ikenm2, ikstp1
75      !! * Local declarations
76      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
77      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
78         zwd, zws, zwi,          &  ! ???
79         zwx, zwy, zwt              ! ???
80      REAL(wp) ::  ztra      ! temporary scalars
81
82      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   &
83         ztrd
84      CHARACTER (len=22) :: charout
85      !!---------------------------------------------------------------------
86
87      IF( kt == nittrc000 ) THEN
88         WRITE(numout,*)
89         WRITE(numout,*) 'trc_zdf_implicit : vertical tracer mixing'
90         WRITE(numout,*) '~~~~~~~~~~~~~~~'
91      ENDIF
92
93      ! 0. Local constant initialization
94      ! --------------------------------
95      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
96         ! time step = 2 rdttra with Arakawa or TVD advection scheme
97         IF( neuler == 0 .AND. kt == nittrc000 ) THEN
98            rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
99         ELSEIF( kt <= nittrc000 + ndttrc ) THEN
100            rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
101         ENDIF
102      ELSE
103         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)     
104      ENDIF
105
106      DO jn = 1 , jptra
107
108    ! Initialisation     
109    zwd( 1 ,:,:) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
110    zws( 1 ,:,:) = 0.e0     ;     zws(jpi,:,:) = 0.e0
111    zwi( 1 ,:,:) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
112    zwt( 1 ,:,:) = 0.e0     ;     zwt(jpi,:,:) = 0.e0     
113         zwt(  :,:,1) = 0.e0     ;     zwt(  :,:,jpk) = 0.e0
114         !                                         
115         ! 0. Matrix construction
116         ! ----------------------
117
118         ! Diagonal, inferior, superior
119         ! (including the bottom boundary condition via avs masked
120         DO jk = 1, jpkm1                                                     
121            DO jj = 2, jpjm1                                     
122               DO ji = fs_2, fs_jpim1   ! vector opt.
123                  zwi(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  ) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
124                  zws(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
125                  zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
126               END DO
127            END DO
128         END DO
129
130         ! Surface boudary conditions
131         DO jj = 2, jpjm1       
132            DO ji = fs_2, fs_jpim1
133               zwi(ji,jj,1) = 0.e0
134               zwd(ji,jj,1) = 1. - zws(ji,jj,1)
135            END DO
136         END DO
137         
138         ! Second member construction
139         DO jk = 1, jpkm1
140            DO jj = 2, jpjm1     
141               DO ji = fs_2, fs_jpim1
142                  zwy(ji,jj,jk) = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)
143               END DO
144            END DO
145         END DO
146         
147         
148         ! Matrix inversion from the first level
149         !----------------------------------------------------------------------
150         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
151         !
152         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
153         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
154         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
155         !        (        ...               )( ...  ) ( ...  )
156         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
157         !
158         !   m is decomposed in the product of an upper and lower triangular matrix
159         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
160         !   The second member is in 2d array zwy
161         !   The solution is in 2d array zwx
162         !   The 3d arry zwt is a work space array
163         !   zwy is used and then used as a work space array : its value is modified!
164         !
165         !   N.B. the starting vertical index (ikst) is equal to 1 except for
166         !   the resolution of tke matrix where surface tke value is prescribed
167         !   so that ikstrt=2.
168        ikst = 1
169        ikstp1 = ikst + 1
170        ikenm2 = jpk - 2
171       
172        DO jj = 2, jpjm1
173           DO ji = fs_2, fs_jpim1
174              zwt(ji,jj,ikst) = zwd(ji,jj,ikst)
175           END DO
176        END DO
177       
178        DO jk = ikstp1, jpkm1
179           DO jj = 2, jpjm1
180              DO ji = fs_2, fs_jpim1
181                 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
182              END DO
183           END DO
184        END DO
185       
186        DO jk = ikstp1, jpkm1
187           DO jj = 2, jpjm1
188              DO ji = fs_2, fs_jpim1
189                 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * zwy(ji,jj,jk-1)
190              END DO
191           END DO
192        END DO
193       
194        DO jj = 2, jpjm1
195           DO ji = fs_2, fs_jpim1
196              zwx(ji,jj,jpkm1) = zwy(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
197           END DO
198        END DO
199       
200        DO jk = ikenm2, ikst, -1
201           DO jj = 2, jpjm1
202              DO ji = fs_2, fs_jpim1
203                 zwx(ji,jj,jk) = ( zwy(ji,jj,jk) - zws(ji,jj,jk) * zwx(ji,jj,jk+1) ) / zwt(ji,jj,jk)
204              END DO
205           END DO
206        END DO
207
208         
209#if defined key_trc_diatrd
210         ! Compute and save the vertical diffusive of tracers trends
211#  if defined key_trcldf_iso
212         DO jk = 1, jpkm1
213            DO jj = 2, jpjm1
214               DO ji = fs_2, fs_jpim1   ! vector opt.
215                  ztra = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
216                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn) + trtrd(ji,jj,jk,ikeep(jn),6)
217               END DO
218            END DO
219         END DO
220#  else
221         DO jk = 1, jpkm1
222            DO jj = 2, jpjm1
223               DO ji = fs_2, fs_jpim1   ! vector opt.
224                  ztra = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
225                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn)
226               END DO
227            END DO
228         END DO
229#  endif
230#endif 
231         ! Save the masked passive tracer after in tra
232         ! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
233         !                  it will not be done in tranxt)
234         DO jk = 1, jpkm1
235            DO jj = 2, jpjm1
236               DO ji = fs_2, fs_jpim1
237                  tra(ji,jj,jk,jn) = zwx(ji,jj,jk) * tmask(ji,jj,jk)
238               END DO
239            END DO
240         END DO
241
242         IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
243            ztrd(:,:,:,:) = 0.
244            DO jk = 1, jpkm1
245               DO jj = 2, jpjm1
246                  DO ji = fs_2, fs_jpim1
247                     ztrd(ji,jj,jk,jn) = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
248                  END DO
249               END DO
250            END DO
251         ENDIF
252
253      END DO
254
255      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
256         WRITE(charout, FMT="('zdf - imp')")
257         CALL prt_ctl_trc_info(charout)
258         CALL prt_ctl_trc(tab4d=ztrd, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
259      ENDIF
260
261   END SUBROUTINE trc_zdf_imp
262
263#else
264   !!----------------------------------------------------------------------
265   !!   Dummy module :                      NO passive tracer
266   !!----------------------------------------------------------------------
267CONTAINS
268   SUBROUTINE trc_zdf_imp (kt )              ! Empty routine
269      INTEGER, INTENT(in) :: kt
270      WRITE(*,*) 'trc_zdf_imp: You should not have seen this print! error?', kt
271   END SUBROUTINE trc_zdf_imp
272#endif
273   
274!!==============================================================================
275END MODULE trczdf_imp
Note: See TracBrowser for help on using the repository browser.