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

Last change on this file since 1438 was 1438, checked in by rblod, 15 years ago

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