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

source: branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2205

Last change on this file since 2205 was 2205, checked in by acc, 14 years ago

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