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 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 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   !!                 Vector optimization, use k-j-i loops.
20   !!----------------------------------------------------------------------
21   !! * Modules used
22   USE oce             ! ocean dynamics and tracers variables
23   USE dom_oce         ! ocean space and time domain variables
24   USE zdf_oce         ! ocean vertical physics variables
25   USE ldftra_oce      ! ocean active tracers: lateral physics
26   USE ldfslp          ! lateral physics: slope of diffusion
27   USE trdmod          ! ocean active tracers trends
28   USE trdmod_oce      ! ocean variables trends
29   USE zdfddm          ! ocean vertical physics: double diffusion
30   USE in_out_manager  ! I/O manager
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE prtctl          ! Print control
33   USE domvvl          ! variable volume
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 (k-j-i loops)'
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      DO jk = 2, jpkm1
142         DO jj = 2, jpjm1
143            DO ji = fs_2, fs_jpim1   ! vector opt.
144               zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
145                  & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
146                  &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
147               zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
148# if defined key_zdfddm
149               zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity
150# endif
151            END DO
152         END DO
153      END DO
154
155#else
156      ! No isopycnal diffusion
157      zwt(:,:,:) = avt(:,:,:)
158# if defined key_zdfddm
159      zavsi(:,:,:) = avs(:,:,:)
160# endif
161
162#endif
163
164      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
165      DO jk = 1, jpkm1
166         DO jj = 2, jpjm1
167            DO ji = fs_2, fs_jpim1   ! vector opt.
168               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
169               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                ! after scale factor at T-point
170               ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                        ! now   scale factor at T-point
171               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
172               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
173               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
174            END DO
175         END DO
176      END DO
177
178      ! Surface boudary conditions
179      DO jj = 2, jpjm1
180         DO ji = fs_2, fs_jpim1   ! vector opt.
181            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
182            ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                   ! after scale factor at T-point
183            zwi(ji,jj,1) = 0.e0
184            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
185         END DO
186      END DO
187
188
189      ! II.1. Vertical diffusion on t
190      ! ---------------------------
191
192      !! Matrix inversion from the first level
193      !!----------------------------------------------------------------------
194      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
195      !
196      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
197      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
198      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
199      !        (        ...               )( ...  ) ( ...  )
200      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
201      !
202      !   m is decomposed in the product of an upper and lower triangular matrix
203      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
204      !   The second member is in 2d array zwy
205      !   The solution is in 2d array zwx
206      !   The 3d arry zwt is a work space array
207      !   zwy is used and then used as a work space array : its value is modified!
208
209      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
210      DO jj = 2, jpjm1
211         DO ji = fs_2, fs_jpim1
212            zwt(ji,jj,1) = zwd(ji,jj,1)
213         END DO
214      END DO
215      DO jk = 2, jpkm1
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1
218               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
219            END DO
220         END DO
221      END DO
222
223      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
224      DO jj = 2, jpjm1
225         DO ji = fs_2, fs_jpim1
226            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
227            ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
228            ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1)
229            ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
230         END DO
231      END DO
232      DO jk = 2, jpkm1
233         DO jj = 2, jpjm1
234            DO ji = fs_2, fs_jpim1
235               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
236               ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
237               ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk)
238               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side
239               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
240            END DO
241         END DO
242      END DO
243
244      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
245      ! Save the masked temperature after in ta
246      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
247      DO jj = 2, jpjm1
248         DO ji = fs_2, fs_jpim1
249            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
250         END DO
251      END DO
252      DO jk = jpk-2, 1, -1
253         DO jj = 2, jpjm1
254            DO ji = fs_2, fs_jpim1
255               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)
256            END DO
257         END DO
258      END DO
259
260      ! II.2 Vertical diffusion on salinity
261      ! -----------------------------------
262
263#if defined key_zdfddm
264      ! Rebuild the Matrix as avt /= avs
265
266      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
267      DO jk = 1, jpkm1
268         DO jj = 2, jpjm1
269            DO ji = fs_2, fs_jpim1   ! vector opt.
270               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
271               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                      ! after scale factor at T-point
272               ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                              ! now   scale factor at T-point
273               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
274               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
275               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
276            END DO
277         END DO
278      END DO
279
280      ! Surface boudary conditions
281      DO jj = 2, jpjm1
282         DO ji = fs_2, fs_jpim1   ! vector opt.
283            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
284            ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                          ! after scale factor at T-point
285            zwi(ji,jj,1) = 0.e0
286            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
287         END DO
288      END DO
289#endif
290
291
292      !! Matrix inversion from the first level
293      !!----------------------------------------------------------------------
294      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
295      !
296      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
297      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
298      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
299      !        (        ...               )( ...  ) ( ...  )
300      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
301      !
302      !   m is decomposed in the product of an upper and lower triangular
303      !   matrix
304      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
305      !   The second member is in 2d array zwy
306      !   The solution is in 2d array zwx
307      !   The 3d arry zwt is a work space array
308      !   zwy is used and then used as a work space array : its value is modified!
309
310      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
311      DO jj = 2, jpjm1
312         DO ji = fs_2, fs_jpim1
313            zwt(ji,jj,1) = zwd(ji,jj,1)
314         END DO
315      END DO
316      DO jk = 2, jpkm1
317         DO jj = 2, jpjm1
318            DO ji = fs_2, fs_jpim1
319               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
320            END DO
321         END DO
322      END DO
323
324      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
325      DO jj = 2, jpjm1
326         DO ji = fs_2, fs_jpim1
327            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
328            ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                               ! before scale factor at T-point
329            ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1)                        ! now    scale factor at T-point
330            sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
331         END DO
332      END DO
333      DO jk = 2, jpkm1
334         DO jj = 2, jpjm1
335            DO ji = fs_2, fs_jpim1
336               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
337               ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl                            ! before scale factor at T-point
338               ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk)                    ! now    scale factor at T-point
339               zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side
340               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
341            END DO
342         END DO
343      END DO
344
345      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
346      ! Save the masked temperature after in ta
347      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
348      DO jj = 2, jpjm1
349         DO ji = fs_2, fs_jpim1
350            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
351         END DO
352      END DO
353      DO jk = jpk-2, 1, -1
354         DO jj = 2, jpjm1
355            DO ji = fs_2, fs_jpim1
356               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)
357            END DO
358         END DO
359      END DO
360
361   END SUBROUTINE tra_zdf_imp
362
363   !!==============================================================================
364END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.