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

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

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