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

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

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