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

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

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