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

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 16.2 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: trazdf_imp.F90 1517 2009-07-21 15:41:35Z ctlod $
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      ELSE                         ! no rotation but key_ldfslp defined
153         zwt  (:,:,:) = avt(:,:,:)
154# if defined key_zdfddm
155         zavsi(:,:,:) = avs(:,:,:)       ! avs /= avt (double diffusion mixing)
156# endif
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               ze3ta =  ( 1. - znvvl )   &                            ! after scale factor at T-point
172                  &   +        znvvl    * fse3t_a(ji,jj,jk) 
173               ze3tn =    znvvl          &                            ! now   scale factor at T-point
174                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
175               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
176               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
177               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
178            END DO
179         END DO
180      END DO
181
182      ! Surface boudary conditions
183      DO jj = 2, jpjm1
184         DO ji = fs_2, fs_jpim1   ! vector opt.
185            ze3ta = ( 1. - znvvl )        &                           ! after scale factor at T-point
186               &   +       znvvl      * fse3t_a(ji,jj,1) 
187            zwi(ji,jj,1) = 0.e0
188            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
189         END DO
190      END DO
191
192
193      ! II.1. Vertical diffusion on t
194      ! ---------------------------
195
196      !! Matrix inversion from the first level
197      !!----------------------------------------------------------------------
198      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
199      !
200      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
201      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
202      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
203      !        (        ...               )( ...  ) ( ...  )
204      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
205      !
206      !   m is decomposed in the product of an upper and lower triangular matrix
207      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
208      !   The second member is in 2d array zwy
209      !   The solution is in 2d array zwx
210      !   The 3d arry zwt is a work space array
211      !   zwy is used and then used as a work space array : its value is modified!
212
213      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
214      DO jj = 2, jpjm1
215         DO ji = fs_2, fs_jpim1
216            zwt(ji,jj,1) = zwd(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               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
223            END DO
224         END DO
225      END DO
226
227      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
228      DO jj = 2, jpjm1
229         DO ji = fs_2, fs_jpim1
230            ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
231            ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
232            ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
233         END DO
234      END DO
235      DO jk = 2, jpkm1
236         DO jj = 2, jpjm1
237            DO ji = fs_2, fs_jpim1
238               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
239               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
240               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side
241               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1)
242            END DO
243         END DO
244      END DO
245
246      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
247      ! Save the masked temperature after in ta
248      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
249      DO jj = 2, jpjm1
250         DO ji = fs_2, fs_jpim1
251            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
252         END DO
253      END DO
254      DO jk = jpk-2, 1, -1
255         DO jj = 2, jpjm1
256            DO ji = fs_2, fs_jpim1
257               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)
258            END DO
259         END DO
260      END DO
261
262      ! II.2 Vertical diffusion on salinity
263      ! -----------------------------------
264
265#if defined key_zdfddm
266      ! Rebuild the Matrix as avt /= avs
267
268      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
269      DO jk = 1, jpkm1
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1   ! vector opt.
272               ze3ta =  ( 1. - znvvl )                        &         ! after scale factor at T-point
273                  &   +        znvvl   * fse3t_a(ji,jj,jk)           
274               ze3tn =    znvvl                               &         ! now   scale factor at T-point
275                  &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)
276               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
277               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
278               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
279            END DO
280         END DO
281      END DO
282
283      ! Surface boudary conditions
284      DO jj = 2, jpjm1
285         DO ji = fs_2, fs_jpim1   ! vector opt.
286            ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1)
287            zwi(ji,jj,1) = 0.e0
288            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
289         END DO
290      END DO
291#endif
292
293
294      !! Matrix inversion from the first level
295      !!----------------------------------------------------------------------
296      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
297      !
298      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
299      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
300      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
301      !        (        ...               )( ...  ) ( ...  )
302      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
303      !
304      !   m is decomposed in the product of an upper and lower triangular
305      !   matrix
306      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
307      !   The second member is in 2d array zwy
308      !   The solution is in 2d array zwx
309      !   The 3d arry zwt is a work space array
310      !   zwy is used and then used as a work space array : its value is modified!
311
312      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
313      DO jj = 2, jpjm1
314         DO ji = fs_2, fs_jpim1
315            zwt(ji,jj,1) = zwd(ji,jj,1)
316         END DO
317      END DO
318      DO jk = 2, jpkm1
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1
321               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
322            END DO
323         END DO
324      END DO
325
326      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
327      DO jj = 2, jpjm1
328         DO ji = fs_2, fs_jpim1
329            ze3tb = ( 1. - znvvl )   &                                           ! before scale factor at T-point
330               &   +  znvvl       * fse3t_b(ji,jj,1)
331            ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,1)                    ! now    scale factor at T-point
332            sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
333         END DO
334      END DO
335      DO jk = 2, jpkm1
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1
338               ze3tb = ( 1. - znvvl )   &                                        ! before scale factor at T-point
339                  &   +  znvvl       * fse3t_b(ji,jj,jk)
340               ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)                ! now    scale factor at T-point
341               zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side
342               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
343            END DO
344         END DO
345      END DO
346
347      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
348      ! Save the masked temperature after in ta
349      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
350      DO jj = 2, jpjm1
351         DO ji = fs_2, fs_jpim1
352            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
353         END DO
354      END DO
355      DO jk = jpk-2, 1, -1
356         DO jj = 2, jpjm1
357            DO ji = fs_2, fs_jpim1
358               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)
359            END DO
360         END DO
361      END DO
362      !
363   END SUBROUTINE tra_zdf_imp
364
365   !!==============================================================================
366END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.