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 branches/TAM_V3_0/NEMO/OPA_SRC/TRA – NEMO

source: branches/TAM_V3_0/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 1884

Last change on this file since 1884 was 1884, checked in by rblod, 14 years ago

Light adaptation of NEMO direct model routine to handle TAM

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