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.
trazdf_imp.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 902

Last change on this file since 902 was 902, checked in by rblod, 16 years ago

Correct a bug with key_ldf_slp and ln_traldf_iso = .FALSE., see ticket #118
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.2 KB
Line 
1MODULE trazdf_imp
2   !!==============================================================================
3   !!                 ***  MODULE  trazdf_imp  ***
4   !! Ocean active tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6   !! History :
7   !!   6.0  !  90-10  (B. Blanke)  Original code
8   !!   7.0  !  91-11 (G. Madec)
9   !!        !  92-06 (M. Imbard) correction on tracer trend loops
10   !!        !  96-01 (G. Madec) statement function for e3
11   !!        !  97-05 (G. Madec) vertical component of isopycnal
12   !!        !  97-07 (G. Madec) geopotential diffusion in s-coord
13   !!        !  00-08 (G. Madec) double diffusive mixing
14   !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
15   !!   9.0  !  06-11 (G. Madec)  New step reorganisation
16   !!----------------------------------------------------------------------
17   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
18   !!                 part of the mixing tensor.
19   !!----------------------------------------------------------------------
20   !! * Modules used
21   USE oce             ! ocean dynamics and tracers variables
22   USE dom_oce         ! ocean space and time domain variables
23   USE zdf_oce         ! ocean vertical physics variables
24   USE ldftra_oce      ! ocean active tracers: lateral physics
25   USE ldfslp          ! lateral physics: slope of diffusion
26   USE trdmod          ! ocean active tracers trends
27   USE trdmod_oce      ! ocean variables trends
28   USE zdfddm          ! ocean vertical physics: double diffusion
29   USE in_out_manager  ! I/O manager
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE prtctl          ! Print control
32   USE domvvl          ! variable volume
33   USE traldf          ! lateral mixing type
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC tra_zdf_imp   !  routine called by step.F90
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "ldftra_substitute.h90"
44#  include "zdfddm_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!----------------------------------------------------------------------
48   !!  OPA 9.0 , LOCEAN-IPSL (2005)
49   !!----------------------------------------------------------------------
50CONTAINS
51   
52   SUBROUTINE tra_zdf_imp( kt, p2dt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_zdf_imp  ***
55      !!
56      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
57      !!     including the vertical component of lateral mixing (only for 2nd
58      !!     order operator, for fourth order it is already computed and add
59      !!     to the general trend in traldf.F) and add it to the general trend
60      !!     of the tracer equations.
61      !!
62      !! ** Method  :   The vertical component of the lateral diffusive trends
63      !!      is provided by a 2nd order operator rotated along neutral or geo-
64      !!      potential surfaces to which an eddy induced advection can be
65      !!      added. It is computed using before fields (forward in time) and
66      !!      isopycnal or geopotential slopes computed in routine ldfslp.
67      !!
68      !!      Second part: vertical trend associated with the vertical physics
69      !!      ===========  (including the vertical flux proportional to dk[t]
70      !!                  associated with the lateral mixing, through the
71      !!                  update of avt)
72      !!      The vertical diffusion of tracers (t & s) is given by:
73      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
74      !!      It is computed using a backward time scheme (t=ta).
75      !!      Surface and bottom boundary conditions: no diffusive flux on
76      !!      both tracers (bottom, applied through the masked field avt).
77      !!      Add this trend to the general trend ta,sa :
78      !!         ta = ta + dz( avt dz(t) )
79      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
80      !!
81      !!      Third part: recover avt resulting from the vertical physics
82      !!      ==========  alone, for further diagnostics (for example to
83      !!                  compute the turbocline depth in zdfmxl.F90).
84      !!         avt = zavt
85      !!         (avs = zavs if lk_zdfddm=T )
86      !!
87      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
88      !!
89      !!---------------------------------------------------------------------
90      !! * Modules used
91      USE oce    , ONLY :   zwd   => ua,  &  ! ua used as workspace
92                            zws   => va      ! va   "          "
93      !! * Arguments
94      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
95      REAL(wp), DIMENSION(jpk), INTENT( in ) ::   &
96         p2dt                                ! vertical profile of tracer time-step
97
98      !! * Local declarations
99      INTEGER  ::   ji, jj, jk               ! dummy loop indices
100      REAL(wp) ::   zavi, zrhs, znvvl,     & ! temporary scalars
101         ze3tb, ze3tn, ze3ta, zvsfvvl        ! variable vertical scale factors
102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
103         zwi, zwt, zavsi                     ! workspace arrays
104      !!---------------------------------------------------------------------
105
106      IF( kt == nit000 ) THEN
107         IF(lwp)WRITE(numout,*)
108         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
109         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
110         zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
111      ENDIF
112
113      ! I. Local initialization
114      ! -----------------------
115      zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0
116      zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0
117      zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0
118      zwt  (1,:, : ) = 0.e0     ;     zwt  (jpi,:,:) = 0.e0
119      zavsi(1,:, : ) = 0.e0     ;     zavsi(jpi,:,:) = 0.e0
120      zwt  (:,:,jpk) = 0.e0     ;     zwt  ( : ,:,1) = 0.e0
121      zavsi(:,:,jpk) = 0.e0     ;     zavsi( : ,:,1) = 0.e0
122
123      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
124      ! -------------------
125      IF( lk_vvl ) THEN   ;    znvvl = 1.
126      ELSE                ;    znvvl = 0.e0
127      ENDIF
128
129      ! II. Vertical trend associated with the vertical physics
130      ! =======================================================
131      !     (including the vertical flux proportional to dk[t] associated
132      !      with the lateral mixing, through the avt update)
133      !     dk[ avt dk[ (t,s) ] ] diffusive trends
134
135
136      ! II.0 Matrix construction
137      ! ------------------------
138
139#if defined key_ldfslp
140      ! update and save of avt (and avs if double diffusive mixing)
141      IF( nldf == 1) THEN
142         DO jk = 2, jpkm1
143            DO jj = 2, jpjm1
144               DO ji = fs_2, fs_jpim1   ! vector opt.
145                  zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
146                     & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
147                     &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
148                  zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
149# if defined key_zdfddm
150                  zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity
151# endif
152               END DO
153            END DO
154         END DO
155      ENDIF   
156
157#else
158      ! No isopycnal diffusion
159      zwt(:,:,:) = avt(:,:,:)
160# if defined key_zdfddm
161      zavsi(:,:,:) = avs(:,:,:)
162# endif
163
164#endif
165
166      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
167      DO jk = 1, jpkm1
168         DO jj = 2, jpjm1
169            DO ji = fs_2, fs_jpim1   ! vector opt.
170               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
171               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                ! after scale factor at T-point
172               ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                        ! now   scale factor at T-point
173               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
174               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
175               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
176            END DO
177         END DO
178      END DO
179
180      ! Surface boudary conditions
181      DO jj = 2, jpjm1
182         DO ji = fs_2, fs_jpim1   ! vector opt.
183            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
184            ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                   ! after scale factor at T-point
185            zwi(ji,jj,1) = 0.e0
186            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
187         END DO
188      END DO
189
190
191      ! II.1. Vertical diffusion on t
192      ! ---------------------------
193
194      !! Matrix inversion from the first level
195      !!----------------------------------------------------------------------
196      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
197      !
198      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
199      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
200      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
201      !        (        ...               )( ...  ) ( ...  )
202      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
203      !
204      !   m is decomposed in the product of an upper and lower triangular matrix
205      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
206      !   The second member is in 2d array zwy
207      !   The solution is in 2d array zwx
208      !   The 3d arry zwt is a work space array
209      !   zwy is used and then used as a work space array : its value is modified!
210
211      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
212      DO jj = 2, jpjm1
213         DO ji = fs_2, fs_jpim1
214            zwt(ji,jj,1) = zwd(ji,jj,1)
215         END DO
216      END DO
217      DO jk = 2, jpkm1
218         DO jj = 2, jpjm1
219            DO ji = fs_2, fs_jpim1
220               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
221            END DO
222         END DO
223      END DO
224
225      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
226      DO jj = 2, jpjm1
227         DO ji = fs_2, fs_jpim1
228            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
229            ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
230            ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1)
231            ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
232         END DO
233      END DO
234      DO jk = 2, jpkm1
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1
237               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
238               ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
239               ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk)
240               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side
241               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
242            END DO
243         END DO
244      END DO
245
246      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
247      ! Save the masked temperature after in ta
248      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
249      DO jj = 2, jpjm1
250         DO ji = fs_2, fs_jpim1
251            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
252         END DO
253      END DO
254      DO jk = jpk-2, 1, -1
255         DO jj = 2, jpjm1
256            DO ji = fs_2, fs_jpim1
257               ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
258            END DO
259         END DO
260      END DO
261
262      ! II.2 Vertical diffusion on salinity
263      ! -----------------------------------
264
265#if defined key_zdfddm
266      ! Rebuild the Matrix as avt /= avs
267
268      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
269      DO jk = 1, jpkm1
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1   ! vector opt.
272               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
273               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                      ! after scale factor at T-point
274               ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                              ! now   scale factor at T-point
275               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
276               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
277               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
278            END DO
279         END DO
280      END DO
281
282      ! Surface boudary conditions
283      DO jj = 2, jpjm1
284         DO ji = fs_2, fs_jpim1   ! vector opt.
285            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
286            ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                          ! after scale factor at T-point
287            zwi(ji,jj,1) = 0.e0
288            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
289         END DO
290      END DO
291#endif
292
293
294      !! Matrix inversion from the first level
295      !!----------------------------------------------------------------------
296      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
297      !
298      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
299      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
300      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
301      !        (        ...               )( ...  ) ( ...  )
302      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
303      !
304      !   m is decomposed in the product of an upper and lower triangular
305      !   matrix
306      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
307      !   The second member is in 2d array zwy
308      !   The solution is in 2d array zwx
309      !   The 3d arry zwt is a work space array
310      !   zwy is used and then used as a work space array : its value is modified!
311
312      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
313      DO jj = 2, jpjm1
314         DO ji = fs_2, fs_jpim1
315            zwt(ji,jj,1) = zwd(ji,jj,1)
316         END DO
317      END DO
318      DO jk = 2, jpkm1
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1
321               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
322            END DO
323         END DO
324      END DO
325
326      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
327      DO jj = 2, jpjm1
328         DO ji = fs_2, fs_jpim1
329            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
330            ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                               ! before scale factor at T-point
331            ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1)                        ! now    scale factor at T-point
332            sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
333         END DO
334      END DO
335      DO jk = 2, jpkm1
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1
338               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
339               ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                            ! before scale factor at T-point
340               ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk)                    ! now    scale factor at T-point
341               zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side
342               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
343            END DO
344         END DO
345      END DO
346
347      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
348      ! Save the masked temperature after in ta
349      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
350      DO jj = 2, jpjm1
351         DO ji = fs_2, fs_jpim1
352            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
353         END DO
354      END DO
355      DO jk = jpk-2, 1, -1
356         DO jj = 2, jpjm1
357            DO ji = fs_2, fs_jpim1
358               sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
359            END DO
360         END DO
361      END DO
362
363   END SUBROUTINE tra_zdf_imp
364
365   !!==============================================================================
366END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.