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

Last change on this file since 527 was 503, checked in by opalod, 18 years ago

nemo_v1_update_064 : CT : general trends update including the addition of mean windows analysis possibility in the mixed layer
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.6 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
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               ! temporary scalars
100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
101         zwi, zwt, zavsi                     ! workspace arrays
102      !!---------------------------------------------------------------------
103
104      IF( kt == nit000 ) THEN
105         IF(lwp)WRITE(numout,*)
106         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing (k-j-i loops)'
107         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
108         zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
109      ENDIF
110
111      ! I. Local initialization
112      ! -----------------------
113      zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0
114      zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0
115      zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0
116      zwt  (1,:, : ) = 0.e0     ;     zwt  (jpi,:,:) = 0.e0
117      zavsi(1,:, : ) = 0.e0     ;     zavsi(jpi,:,:) = 0.e0
118      zwt  (:,:,jpk) = 0.e0     ;     zwt  ( : ,:,1) = 0.e0
119      zavsi(:,:,jpk) = 0.e0     ;     zavsi( : ,:,1) = 0.e0
120
121      ! II. Vertical trend associated with the vertical physics
122      ! =======================================================
123      !     (including the vertical flux proportional to dk[t] associated
124      !      with the lateral mixing, through the avt update)
125      !     dk[ avt dk[ (t,s) ] ] diffusive trends
126
127
128      ! II.0 Matrix construction
129      ! ------------------------
130
131#if defined key_ldfslp
132      ! update and save of avt (and avs if double diffusive mixing)
133      DO jk = 2, jpkm1
134         DO jj = 2, jpjm1
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
137                  & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
138                  &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
139               zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
140# if defined key_zdfddm
141               zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity
142# endif
143            END DO
144         END DO
145      END DO
146
147      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
148      DO jk = 1, jpkm1
149         DO jj = 2, jpjm1
150            DO ji = fs_2, fs_jpim1   ! vector opt.
151               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
152               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
153               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
154            END DO
155         END DO
156      END DO
157#else
158      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
159      DO jk = 1, jpkm1
160         DO jj = 2, jpjm1
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zwi(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
163               zws(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
164               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
165            END DO
166         END DO
167      END DO
168#endif
169
170      ! Surface boudary conditions
171      DO jj = 2, jpjm1
172         DO ji = fs_2, fs_jpim1   ! vector opt.
173            zwi(ji,jj,1) = 0.e0
174            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
175         END DO
176      END DO
177
178
179      ! II.1. Vertical diffusion on t
180      ! ---------------------------
181
182      !! Matrix inversion from the first level
183      !!----------------------------------------------------------------------
184      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
185      !
186      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
187      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
188      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
189      !        (        ...               )( ...  ) ( ...  )
190      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
191      !
192      !   m is decomposed in the product of an upper and lower triangular matrix
193      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
194      !   The second member is in 2d array zwy
195      !   The solution is in 2d array zwx
196      !   The 3d arry zwt is a work space array
197      !   zwy is used and then used as a work space array : its value is modified!
198
199      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
200      DO jj = 2, jpjm1
201         DO ji = fs_2, fs_jpim1
202            zwt(ji,jj,1) = zwd(ji,jj,1)
203         END DO
204      END DO
205      DO jk = 2, jpkm1
206         DO jj = 2, jpjm1
207            DO ji = fs_2, fs_jpim1
208               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
209            END DO
210         END DO
211      END DO
212
213      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
214      DO jj = 2, jpjm1
215         DO ji = fs_2, fs_jpim1
216            ta(ji,jj,1) = tb(ji,jj,1) + p2dt(1) * ta(ji,jj,1)
217         END DO
218      END DO
219      DO jk = 2, jpkm1
220         DO jj = 2, jpjm1
221            DO ji = fs_2, fs_jpim1
222               zrhs = tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)   ! zrhs=right hand side
223               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
224            END DO
225         END DO
226      END DO
227
228      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
229      ! Save the masked temperature after in ta
230      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
231      DO jj = 2, jpjm1
232         DO ji = fs_2, fs_jpim1
233            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
234         END DO
235      END DO
236      DO jk = jpk-2, 1, -1
237         DO jj = 2, jpjm1
238            DO ji = fs_2, fs_jpim1
239               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)
240            END DO
241         END DO
242      END DO
243
244      ! II.2 Vertical diffusion on salinity
245      ! -----------------------------------
246
247#if defined key_zdfddm
248      ! Rebuild the Matrix as avt /= avs
249
250      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
251# if defined key_ldfslp
252      DO jk = 1, jpkm1
253         DO jj = 2, jpjm1
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
256               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
257               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
258            END DO
259         END DO
260      END DO
261# else
262      DO jk = 1, jpkm1
263         DO jj = 2, jpjm1
264            DO ji = fs_2, fs_jpim1   ! vector opt.
265               zwi(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
266               zws(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
267               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
268            END DO
269         END DO
270      END DO
271# endif
272
273      ! Surface boudary conditions
274      DO jj = 2, jpjm1
275         DO ji = fs_2, fs_jpim1   ! vector opt.
276            zwi(ji,jj,1) = 0.e0
277            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
278         END DO
279      END DO
280#endif
281
282
283      !! Matrix inversion from the first level
284      !!----------------------------------------------------------------------
285      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
286      !
287      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
288      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
289      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
290      !        (        ...               )( ...  ) ( ...  )
291      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
292      !
293      !   m is decomposed in the product of an upper and lower triangular
294      !   matrix
295      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
296      !   The second member is in 2d array zwy
297      !   The solution is in 2d array zwx
298      !   The 3d arry zwt is a work space array
299      !   zwy is used and then used as a work space array : its value is modified!
300
301      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
302      DO jj = 2, jpjm1
303         DO ji = fs_2, fs_jpim1
304            zwt(ji,jj,1) = zwd(ji,jj,1)
305         END DO
306      END DO
307      DO jk = 2, jpkm1
308         DO jj = 2, jpjm1
309            DO ji = fs_2, fs_jpim1
310               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
311            END DO
312         END DO
313      END DO
314
315      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
316      DO jj = 2, jpjm1
317         DO ji = fs_2, fs_jpim1
318            sa(ji,jj,1) = sb(ji,jj,1) + p2dt(1) * sa(ji,jj,1)
319         END DO
320      END DO
321      DO jk = 2, jpkm1
322         DO jj = 2, jpjm1
323            DO ji = fs_2, fs_jpim1
324               zrhs = sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)   ! zrhs=right hand side
325               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
326            END DO
327         END DO
328      END DO
329
330      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
331      ! Save the masked temperature after in ta
332      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
333      DO jj = 2, jpjm1
334         DO ji = fs_2, fs_jpim1
335            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
336         END DO
337      END DO
338      DO jk = jpk-2, 1, -1
339         DO jj = 2, jpjm1
340            DO ji = fs_2, fs_jpim1
341               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)
342            END DO
343         END DO
344      END DO
345
346   END SUBROUTINE tra_zdf_imp
347
348   !!==============================================================================
349END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.