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

Last change on this file since 467 was 457, checked in by opalod, 18 years ago

nemo_v1_update_049:RB: reorganization of tracers part, remove traadv_cen2_atsk.h90 traldf_iso_zps.F90 trazdf_iso.F90 trazdf_iso_vopt.F90, change atsk routines to jki
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.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 :
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  !  04-08 (C. Talandier) New trends organization
16   !!   " "  !  06-11 (G. Madec)  New step reorganisation
17   !!----------------------------------------------------------------------
18   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
19   !!                 part of the mixing tensor.
20   !!                 Vector optimization, use k-j-i loops.
21   !!----------------------------------------------------------------------
22   !! * Modules used
23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain variables
25   USE zdf_oce         ! ocean vertical physics variables
26   USE ldftra_oce      ! ocean active tracers: lateral physics
27   USE ldfslp          ! lateral physics: slope of diffusion
28   USE trdmod          ! ocean active tracers trends
29   USE trdmod_oce      ! ocean variables trends
30   USE zdfddm          ! ocean vertical physics: double diffusion
31   USE in_out_manager  ! I/O manager
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE prtctl          ! Print control
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   !!----------------------------------------------------------------------
50CONTAINS
51   
52   SUBROUTINE tra_zdf_imp( kt, p2dt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_zdf_imp  ***
55      !!
56      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
57      !!     including the vertical component of lateral mixing (only for 2nd
58      !!     order operator, for fourth order it is already computed and add
59      !!     to the general trend in traldf.F) and add it to the general trend
60      !!     of the tracer equations.
61      !!
62      !! ** Method  :   The vertical component of the lateral diffusive trends
63      !!      is provided by a 2nd order operator rotated along neural or geo-
64      !!      potential surfaces to which an eddy induced advection can be
65      !!      added. It is computed using before fields (forward in time) and
66      !!      isopycnal or geopotential slopes computed in routine ldfslp.
67      !!
68      !!      Second part: vertical trend associated with the vertical physics
69      !!      ===========  (including the vertical flux proportional to dk[t]
70      !!                  associated with the lateral mixing, through the
71      !!                  update of avt)
72      !!      The vertical diffusion of tracers (t & s) is given by:
73      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
74      !!      It is computed using a backward time scheme (t=ta).
75      !!      Surface and bottom boundary conditions: no diffusive flux on
76      !!      both tracers (bottom, applied through the masked field avt).
77      !!      Add this trend to the general trend ta,sa :
78      !!         ta = ta + dz( avt dz(t) )
79      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
80      !!
81      !!      Third part: recover avt resulting from the vertical physics
82      !!      ==========  alone, for further diagnostics (for example to
83      !!                  compute the turbocline depth in zdfmxl.F90).
84      !!         avt = zavt
85      !!         (avs = zavs if lk_zdfddm=T )
86      !!
87      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
88      !!
89      !!---------------------------------------------------------------------
90      !! * Modules used
91      USE oce    , ONLY :   zwd   => ua,  &  ! ua used as workspace
92                            zws   => va      ! va   "          "
93      !! * Arguments
94      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
95      REAL(wp), DIMENSION(jpk), INTENT( in ) ::   &
96         p2dt                                ! vertical profile of tracer time-step
97
98      !! * Local declarations
99      INTEGER  ::   ji, jj, jk               ! dummy loop indices
100      REAL(wp) ::   zavi, zrhs               ! temporary scalars
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 (k-j-i loops)'
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      ! II. Vertical trend associated with the vertical physics
123      ! =======================================================
124      !     (including the vertical flux proportional to dk[t] associated
125      !      with the lateral mixing, through the avt update)
126      !     dk[ avt dk[ (t,s) ] ] diffusive trends
127
128
129      ! II.0 Matrix construction
130      ! ------------------------
131
132#if defined key_ldfslp
133      ! update and save of avt (and avs if double diffusive mixing)
134      DO jk = 2, jpkm1
135         DO jj = 2, jpjm1
136            DO ji = fs_2, fs_jpim1   ! vector opt.
137               zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
138                  & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
139                  &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
140               zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
141# if defined key_zdfddm
142               zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity
143# endif
144            END DO
145         END DO
146      END DO
147
148      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
149      DO jk = 1, jpkm1
150         DO jj = 2, jpjm1
151            DO ji = fs_2, fs_jpim1   ! vector opt.
152               zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
153               zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
154               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
155            END DO
156         END DO
157      END DO
158#else
159      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
160      DO jk = 1, jpkm1
161         DO jj = 2, jpjm1
162            DO ji = fs_2, fs_jpim1   ! vector opt.
163               zwi(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
164               zws(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
165               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
166            END DO
167         END DO
168      END DO
169#endif
170
171      ! Surface boudary conditions
172      DO jj = 2, jpjm1
173         DO ji = fs_2, fs_jpim1   ! vector opt.
174            zwi(ji,jj,1) = 0.e0
175            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
176         END DO
177      END DO
178
179
180      ! II.1. Vertical diffusion on t
181      ! ---------------------------
182
183      !! Matrix inversion from the first level
184      !!----------------------------------------------------------------------
185      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
186      !
187      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
188      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
189      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
190      !        (        ...               )( ...  ) ( ...  )
191      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
192      !
193      !   m is decomposed in the product of an upper and lower triangular matrix
194      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
195      !   The second member is in 2d array zwy
196      !   The solution is in 2d array zwx
197      !   The 3d arry zwt is a work space array
198      !   zwy is used and then used as a work space array : its value is modified!
199
200      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
201      DO jj = 2, jpjm1
202         DO ji = fs_2, fs_jpim1
203            zwt(ji,jj,1) = zwd(ji,jj,1)
204         END DO
205      END DO
206      DO jk = 2, jpkm1
207         DO jj = 2, jpjm1
208            DO ji = fs_2, fs_jpim1
209               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
210            END DO
211         END DO
212      END DO
213
214      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
215      DO jj = 2, jpjm1
216         DO ji = fs_2, fs_jpim1
217            ta(ji,jj,1) = tb(ji,jj,1) + p2dt(1) * ta(ji,jj,1)
218         END DO
219      END DO
220      DO jk = 2, jpkm1
221         DO jj = 2, jpjm1
222            DO ji = fs_2, fs_jpim1
223               zrhs = tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)   ! zrhs=right hand side
224               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
225            END DO
226         END DO
227      END DO
228
229      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
230      ! Save the masked temperature after in ta
231      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
232      DO jj = 2, jpjm1
233         DO ji = fs_2, fs_jpim1
234            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
235         END DO
236      END DO
237      DO jk = jpk-2, 1, -1
238         DO jj = 2, jpjm1
239            DO ji = fs_2, fs_jpim1
240               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)
241            END DO
242         END DO
243      END DO
244
245      ! II.2 Vertical diffusion on salinity
246      ! -----------------------------------
247
248#if defined key_zdfddm
249      ! Rebuild the Matrix as avt /= avs
250
251      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
252# if defined key_ldfslp
253      DO jk = 1, jpkm1
254         DO jj = 2, jpjm1
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
257               zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
258               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
259            END DO
260         END DO
261      END DO
262# else
263      DO jk = 1, jpkm1
264         DO jj = 2, jpjm1
265            DO ji = fs_2, fs_jpim1   ! vector opt.
266               zwi(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
267               zws(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
268               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
269            END DO
270         END DO
271      END DO
272# endif
273
274      ! Surface boudary conditions
275      DO jj = 2, jpjm1
276         DO ji = fs_2, fs_jpim1   ! vector opt.
277            zwi(ji,jj,1) = 0.e0
278            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
279         END DO
280      END DO
281#endif
282
283
284      !! Matrix inversion from the first level
285      !!----------------------------------------------------------------------
286      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
287      !
288      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
289      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
290      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
291      !        (        ...               )( ...  ) ( ...  )
292      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
293      !
294      !   m is decomposed in the product of an upper and lower triangular
295      !   matrix
296      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
297      !   The second member is in 2d array zwy
298      !   The solution is in 2d array zwx
299      !   The 3d arry zwt is a work space array
300      !   zwy is used and then used as a work space array : its value is modified!
301
302      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
303      DO jj = 2, jpjm1
304         DO ji = fs_2, fs_jpim1
305            zwt(ji,jj,1) = zwd(ji,jj,1)
306         END DO
307      END DO
308      DO jk = 2, jpkm1
309         DO jj = 2, jpjm1
310            DO ji = fs_2, fs_jpim1
311               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
312            END DO
313         END DO
314      END DO
315
316      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
317      DO jj = 2, jpjm1
318         DO ji = fs_2, fs_jpim1
319            sa(ji,jj,1) = sb(ji,jj,1) + p2dt(1) * sa(ji,jj,1)
320         END DO
321      END DO
322      DO jk = 2, jpkm1
323         DO jj = 2, jpjm1
324            DO ji = fs_2, fs_jpim1
325               zrhs = sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)   ! zrhs=right hand side
326               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
327            END DO
328         END DO
329      END DO
330
331      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
332      ! Save the masked temperature after in ta
333      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
334      DO jj = 2, jpjm1
335         DO ji = fs_2, fs_jpim1
336            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
337         END DO
338      END DO
339      DO jk = jpk-2, 1, -1
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1
342               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)
343            END DO
344         END DO
345      END DO
346
347   END SUBROUTINE tra_zdf_imp
348
349   !!==============================================================================
350END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.