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 tags/nemo_v3_1_beta/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v3_1_beta/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 4528

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

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