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_iso_vopt.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90 @ 216

Last change on this file since 216 was 216, checked in by opalod, 19 years ago

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.4 KB
Line 
1MODULE trazdf_iso_vopt
2   !!==============================================================================
3   !!                 ***  MODULE  trazdf_iso_vopt  ***
4   !! Ocean active tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6#if defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
9   !!----------------------------------------------------------------------
10   !!   tra_zdf_iso_vopt : Update the tracer trend with the vertical part of
11   !!                  the isopycnal or geopotential s-coord. operator and
12   !!                  the vertical diffusion. vector optimization, use
13   !!                  k-j-i loops.
14   !!   tra_zdf_iso  :
15   !!   tra_zdf_zdf  :
16   !!----------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE zdf_oce         ! ocean vertical physics variables
21   USE ldftra_oce      ! ocean active tracers: lateral physics
22   USE trdmod          ! ocean active tracers trends
23   USE trdmod_oce      ! ocean variables trends
24   USE ldfslp          ! iso-neutral slopes
25   USE zdfddm          ! ocean vertical physics: double diffusion
26   USE in_out_manager  ! I/O manager
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !! * Routine accessibility
33   PUBLIC tra_zdf_iso_vopt   !  routine called by step.F90
34
35   !! * Module variables
36   REAL(wp), DIMENSION(jpk) ::  &
37      r2dt                          ! vertical profile of 2 x time-step
38   REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &
39      tavg, savg,                     & ! workspace arrays
40      tdta, tdsa                        ! workspace arrays
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "ldftra_substitute.h90"
45#  include "ldfeiv_substitute.h90"
46#  include "zdfddm_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49
50CONTAINS
51   
52   SUBROUTINE tra_zdf_iso_vopt( kt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_zdf_iso_vopt  ***
55      !!
56      !! ** Purpose :
57      !! ** Method  :
58      !! ** Action  :
59      !!
60      !! History :
61      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
62      !!   9.0  !  04-08  (C. Talandier) New trends organization
63      !!---------------------------------------------------------------------
64      !! * Arguments
65      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
66      !! * Local variables
67      REAL(wp) :: zta, zsa               ! temporary scalars
68      !!---------------------------------------------------------------------
69      !!  OPA 8.5, LODYC-IPSL (2002)
70      !!---------------------------------------------------------------------
71
72      IF( kt == nit000 ) THEN
73         IF(lwp)WRITE(numout,*)
74         IF(lwp)WRITE(numout,*) 'tra_zdf_iso_vopt : vertical mixing computation'
75         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~~   is  iso-neutral diffusion : implicit vertical time stepping'
76#if defined key_diaeiv 
77         w_eiv(:,:,:) = 0.e0
78#endif
79      ENDIF
80
81      ! initialization step
82      tavg(:,:,:) = 0.e0
83      savg(:,:,:) = 0.e0
84
85      ! I. vertical extra-diagonal part of the rotated tensor
86      ! -----------------------------------------------------
87
88      CALL tra_zdf_iso( kt )
89
90      IF(l_ctl) THEN         ! print mean trends (used for debugging)
91         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
92         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
93         WRITE(numout,*) ' zdf 1- Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
94         t_ctl = zta   ;   s_ctl = zsa
95      ENDIF
96
97      ! II. vertical diffusion (including the vertical diagonal part of the rotated tensor)
98      ! ----------------------
99
100      CALL tra_zdf_zdf( kt )
101
102      IF(l_ctl) THEN         ! print mean trends (used for debugging)
103         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
104         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
105         WRITE(numout,*) ' zdf 2- Ta: ', zta, ' Sa: ', zsa
106      ENDIF
107
108   END SUBROUTINE tra_zdf_iso_vopt
109
110
111   SUBROUTINE tra_zdf_zdf( kt )
112      !!----------------------------------------------------------------------
113      !!                 ***  ROUTINE tra_zdf_zdf  ***
114      !!                   
115      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
116      !!     including the vertical component of lateral mixing (only for 2nd
117      !!     order operator, for fourth order it is already computed and add
118      !!     to the general trend in traldf.F) and add it to the general trend
119      !!     of the tracer equations.
120      !!
121      !! ** Method  :   The vertical component of the lateral diffusive trends
122      !!      is provided by a 2nd order operator rotated along neural or geo-
123      !!      potential surfaces to which an eddy induced advection can be
124      !!      added. It is computed using before fields (forward in time) and
125      !!      isopycnal or geopotential slopes computed in routine ldfslp.
126      !!
127      !!      Second part: vertical trend associated with the vertical physics
128      !!      ===========  (including the vertical flux proportional to dk[t]
129      !!                  associated with the lateral mixing, through the
130      !!                  update of avt)
131      !!      The vertical diffusion of tracers (t & s) is given by:
132      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
133      !!      It is computed using a backward time scheme (t=ta).
134      !!      Surface and bottom boundary conditions: no diffusive flux on
135      !!      both tracers (bottom, applied through the masked field avt).
136      !!      Add this trend to the general trend ta,sa :
137      !!         ta = ta + dz( avt dz(t) )
138      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
139      !!
140      !!      Third part: recover avt resulting from the vertical physics
141      !!      ==========  alone, for further diagnostics (for example to
142      !!                  compute the turbocline depth in zdfmxl.F90).
143      !!         avt = zavt
144      !!         (avs = zavs if lk_zdfddm=T )
145      !!
146      !!      'key_trdtra' defined: trend saved for futher diagnostics.
147      !!
148      !!      macro-tasked on vertical slab (jj-loop)
149      !!
150      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
151      !!              - Save the trend in (ztdta,ztdsa) ('key_trdtra')
152      !!
153      !! History :
154      !!   6.0  !  90-10  (B. Blanke)  Original code
155      !!   7.0  !  91-11 (G. Madec)
156      !!        !  92-06 (M. Imbard) correction on tracer trend loops
157      !!        !  96-01 (G. Madec) statement function for e3
158      !!        !  97-05 (G. Madec) vertical component of isopycnal
159      !!        !  97-07 (G. Madec) geopotential diffusion in s-coord
160      !!        !  00-08 (G. Madec) double diffusive mixing
161      !!   8.5  !  02-08 (G. Madec)  F90: Free form and module
162      !!   9.0  !  04-08 (C. Talandier) New trends organization
163      !!---------------------------------------------------------------------
164      !! * Modules used
165      USE oce    , ONLY :   zwd   => ua,  &  ! ua, va used as
166                            zws   => va      ! workspace
167      !! * Arguments
168      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
169
170      !! * Local declarations
171      INTEGER ::   ji, jj, jk                ! dummy loop indices
172      REAL(wp) ::   &
173         zavi, zrhs                          ! temporary scalars
174      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
175         zwi, zwt, zavsi                     ! temporary workspace arrays
176      !!---------------------------------------------------------------------
177      !!  OPA 8.5, LODYC-IPSL (2002)
178      !!---------------------------------------------------------------------
179
180
181      ! I. Local constant initialization
182      ! --------------------------------
183      ! ... time step = 2 rdttra ex
184      IF( neuler == 0 .AND. kt == nit000 ) THEN
185         r2dt(:) =  rdttra(:)              ! restarting with Euler time stepping
186      ELSEIF( kt <= nit000 + 1) THEN
187         r2dt(:) = 2. * rdttra(:)          ! leapfrog
188      ENDIF
189
190      ! Save ta and sa trends
191      IF( l_trdtra )   THEN
192         tdta(:,:,:) = ta(:,:,:) 
193         tdsa(:,:,:) = sa(:,:,:) 
194      ENDIF
195
196      zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0
197      zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0
198      zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0
199      zwt  (1,:, : ) = 0.e0     ;     zwt  (jpi,:,:) = 0.e0
200      zavsi(1,:, : ) = 0.e0     ;     zavsi(jpi,:,:) = 0.e0
201      zwt  (:,:,jpk) = 0.e0     ;     zwt  ( : ,:,1) = 0.e0
202      zavsi(:,:,jpk) = 0.e0     ;     zavsi( : ,:,1) = 0.e0
203
204      ! II. Vertical trend associated with the vertical physics
205      ! =======================================================
206      !     (including the vertical flux proportional to dk[t] associated
207      !      with the lateral mixing, through the avt update)
208      !     dk[ avt dk[ (t,s) ] ] diffusive trends
209
210
211      ! II.0 Matrix construction
212      ! ------------------------
213
214      ! update and save of avt (and avs if double diffusive mixing)
215      DO jk = 2, jpkm1
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218               zavi = fsahtw(ji,jj,jk) * (                 &   ! vertical mixing coef. due to lateral mixing
219                    wslpi(ji,jj,jk) * wslpi(ji,jj,jk)      &
220                  + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
221               zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi            ! zwt=avt+zavi (total vertical mixing coef. on temperature)
222#if defined key_zdfddm
223               zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi        ! dd mixing: zavsi = total vertical mixing coef. on salinity
224#endif
225            END DO
226         END DO
227      END DO
228
229      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
230      DO jk = 1, jpkm1
231         DO jj = 2, jpjm1
232            DO ji = fs_2, fs_jpim1   ! vector opt.
233               zwi(ji,jj,jk) = - r2dt(jk) * zwt(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
234               zws(ji,jj,jk) = - r2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
235               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
236            END DO
237         END DO
238      END DO
239
240      ! Surface boudary conditions
241      DO jj = 2, jpjm1
242         DO ji = fs_2, fs_jpim1   ! vector opt.
243            zwi(ji,jj,1) = 0.e0
244            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
245         END DO
246      END DO
247
248
249      ! II.1. Vertical diffusion on t
250      ! ---------------------------
251
252      !! Matrix inversion from the first level
253      !!----------------------------------------------------------------------
254      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
255      !
256      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
257      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
258      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
259      !        (        ...               )( ...  ) ( ...  )
260      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
261      !
262      !   m is decomposed in the product of an upper and lower triangular matrix
263      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
264      !   The second member is in 2d array zwy
265      !   The solution is in 2d array zwx
266      !   The 3d arry zwt is a work space array
267      !   zwy is used and then used as a work space array : its value is modified!
268
269      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
270      DO jj = 2, jpjm1
271         DO ji = fs_2, fs_jpim1
272            zwt(ji,jj,1) = zwd(ji,jj,1)
273         END DO
274      END DO
275      DO jk = 2, jpkm1
276         DO jj = 2, jpjm1
277            DO ji = fs_2, fs_jpim1
278               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
279            END DO
280         END DO
281      END DO
282
283      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
284      DO jj = 2, jpjm1
285         DO ji = fs_2, fs_jpim1
286            ta(ji,jj,1) = tb(ji,jj,1) + r2dt(1) * ta(ji,jj,1)
287         END DO
288      END DO
289      DO jk = 2, jpkm1
290         DO jj = 2, jpjm1
291            DO ji = fs_2, fs_jpim1
292               zrhs = tb(ji,jj,jk) + r2dt(jk) * ta(ji,jj,jk)   ! zrhs=right hand side
293               ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
294            END DO
295         END DO
296      END DO
297
298      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
299      ! Save the masked temperature after in ta
300      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
301      DO jj = 2, jpjm1
302         DO ji = fs_2, fs_jpim1
303            ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
304         END DO
305      END DO
306      DO jk = jpk-2, 1, -1
307         DO jj = 2, jpjm1
308            DO ji = fs_2, fs_jpim1
309               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)
310            END DO
311         END DO
312      END DO
313
314      ! II.2 Vertical diffusion on salinity
315      ! -----------------------------------
316
317#if defined key_zdfddm
318      ! Rebuild the Matrix as avt /= avs
319
320      ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
321      DO jk = 1, jpkm1
322         DO jj = 2, jpjm1
323            DO ji = fs_2, fs_jpim1   ! vector opt.
324               zwi(ji,jj,jk) = - r2dt(jk) * zavsi(ji,jj,jk  ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
325               zws(ji,jj,jk) = - r2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
326               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
327            END DO
328         END DO
329      END DO
330
331      ! Surface boudary conditions
332      DO jj = 2, jpjm1
333         DO ji = fs_2, fs_jpim1   ! vector opt.
334            zwi(ji,jj,1) = 0.e0
335            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
336         END DO
337      END DO
338#endif
339
340
341      !! Matrix inversion from the first level
342      !!----------------------------------------------------------------------
343      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
344      !
345      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
346      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
347      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
348      !        (        ...               )( ...  ) ( ...  )
349      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
350      !
351      !   m is decomposed in the product of an upper and lower triangular
352      !   matrix
353      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
354      !   The second member is in 2d array zwy
355      !   The solution is in 2d array zwx
356      !   The 3d arry zwt is a work space array
357      !   zwy is used and then used as a work space array : its value is modified!
358
359      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
360      DO jj = 2, jpjm1
361         DO ji = fs_2, fs_jpim1
362            zwt(ji,jj,1) = zwd(ji,jj,1)
363         END DO
364      END DO
365      DO jk = 2, jpkm1
366         DO jj = 2, jpjm1
367            DO ji = fs_2, fs_jpim1
368               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
369            END DO
370         END DO
371      END DO
372
373      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
374      DO jj = 2, jpjm1
375         DO ji = fs_2, fs_jpim1
376            sa(ji,jj,1) = sb(ji,jj,1) + r2dt(1) * sa(ji,jj,1)
377         END DO
378      END DO
379      DO jk = 2, jpkm1
380         DO jj = 2, jpjm1
381            DO ji = fs_2, fs_jpim1
382               zrhs = sb(ji,jj,jk) + r2dt(jk) * sa(ji,jj,jk)   ! zrhs=right hand side
383               sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
384            END DO
385         END DO
386      END DO
387
388      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
389      ! Save the masked temperature after in ta
390      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
391      DO jj = 2, jpjm1
392         DO ji = fs_2, fs_jpim1
393            sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
394         END DO
395      END DO
396      DO jk = jpk-2, 1, -1
397         DO jj = 2, jpjm1
398            DO ji = fs_2, fs_jpim1
399               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)
400            END DO
401         END DO
402      END DO
403
404
405      ! save the trends for diagnostic
406      ! compute the vertical diffusive trends in substracting the previous
407      ! trends tdta()/tdsa() to the new one computed via dT/dt or dS/dt
408      ! i.e. with the new temperature/salinity ta/sa computed above
409      IF( l_trdtra )   THEN
410         IF( l_traldf_iso)   THEN
411            DO jk = 1, jpkm1
412               tdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - tdta(:,:,jk) + tavg(:,:,jk) 
413               tdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - tdsa(:,:,jk) + savg(:,:,jk) 
414            END DO
415         ELSE
416            DO jk = 1, jpkm1
417               tdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - tdta(:,:,jk)                             
418               tdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - tdsa(:,:,jk)                             
419            END DO
420         ENDIF
421
422         CALL trd_mod(tdta, tdsa, jpttdzdf, 'TRA', kt)
423      ENDIF
424
425   END SUBROUTINE tra_zdf_zdf
426
427
428   SUBROUTINE tra_zdf_iso( kt )
429      !!----------------------------------------------------------------------
430      !!                  ***  ROUTINE tra_zdf_iso  ***
431      !!
432      !! ** Purpose :
433      !!     Compute the trend due to the vertical tracer diffusion inclu-
434      !!     ding the vertical component of lateral mixing (only for second
435      !!     order operator, for fourth order it is already computed and
436      !!     add to the general trend in traldf.F) and add it to the general
437      !!     trend of the tracer equations.
438      !!
439      !! ** Method :
440      !!         The vertical component of the lateral diffusive trends is
441      !!      provided by a 2nd order operator rotated along neural or geopo-
442      !!      tential surfaces to which an eddy induced advection can be added
443      !!      It is computed using before fields (forward in time) and isopyc-
444      !!      nal or geopotential slopes computed in routine ldfslp.
445      !!
446      !!      First part: vertical trends associated with the lateral mixing
447      !!      ==========  (excluding the vertical flux proportional to dk[t] )
448      !!      vertical fluxes associated with the rotated lateral mixing:
449      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
450      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
451      !!      save avt coef. resulting from vertical physics alone in zavt:
452      !!         zavt = avt
453      !!      update and save in zavt the vertical eddy viscosity coefficient:
454      !!         avt = avt + wslpi^2+wslj^2
455      !!      add vertical Eddy Induced advective fluxes (lk_traldf_eiv=T):
456      !!         zftw = zftw + { di[aht e2u mi(wslpi)]
457      !!                    +dj[aht e1v mj(wslpj)] } mk(tb)
458      !!      take the horizontal divergence of the fluxes:
459      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
460      !!      Add this trend to the general trend (ta,sa):
461      !!         ta = ta + difft
462      !!
463      !! ** Action :
464      !!         Update (ta,sa) arrays with the before vertical diffusion trend
465      !!         Save in (ztdta,ztdsa) arrays the trends if 'key_trdtra' defined
466      !!
467      !! History :
468      !!   6.0  !  90-10  (B. Blanke)  Original code
469      !!   7.0  !  91-11  (G. Madec)
470      !!        !  92-06  (M. Imbard) correction on tracer trend loops
471      !!        !  96-01  (G. Madec) statement function for e3
472      !!        !  97-05  (G. Madec) vertical component of isopycnal
473      !!        !  97-07  (G. Madec) geopotential diffusion in s-coord
474      !!        !  00-08  (G. Madec) double diffusive mixing
475      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
476      !!   9.0  !  04-08  (C. Talandier) New trends organization
477      !!---------------------------------------------------------------------
478      !! * Modules used
479      USE oce    , ONLY :   zwx => ua,  &  ! use ua, va as
480                            zwy => va      ! workspace arrays
481
482      !! * Arguments
483      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
484
485      !! * Local declarations
486      INTEGER ::   ji, jj, jk       ! dummy loop indices
487#if defined key_partial_steps
488      INTEGER ::   iku, ikv
489#endif
490      REAL(wp) ::   &
491         zcoef0, zcoef3,         &  !    "         "
492         zcoef4,                 &  !    "         "
493         zbtr, zmku, zmkv,       &  !    "         "
494#if defined key_traldf_eiv
495         zcoeg3,                 &  !    "         "
496         zuwki, zvwki,           &  !    "         "
497         zuwk, zvwk,             &  !    "         "
498#endif
499         ztav, zsav
500      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
501         zwz, zwt, ztfw, zsfw       ! temporary workspace arrays
502      !!---------------------------------------------------------------------
503      !!  OPA 8.5, LODYC-IPSL (2002)
504      !!---------------------------------------------------------------------
505
506
507      ! 0. Local constant initialization
508      ! --------------------------------
509      zwx (1,:,:) = 0.e0     ;     zwx (jpi,:,:) = 0.e0
510      zwy (1,:,:) = 0.e0     ;     zwy (jpi,:,:) = 0.e0
511      zwz (1,:,:) = 0.e0     ;     zwz (jpi,:,:) = 0.e0
512      zwt (1,:,:) = 0.e0     ;     zwt (jpi,:,:) = 0.e0
513      ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0
514      zsfw(1,:,:) = 0.e0     ;     zsfw(jpi,:,:) = 0.e0
515
516      ! Save ta and sa trends
517      IF( l_trdtra )   THEN
518         tdta(:,:,:) = ta(:,:,:) 
519         tdsa(:,:,:) = sa(:,:,:) 
520      ENDIF
521
522      ! I. Vertical trends associated with lateral mixing
523      ! -------------------------------------------------
524      !    (excluding the vertical flux proportional to dk[t] )
525
526
527      ! I.1 horizontal tracer gradient
528      ! ------------------------------
529
530      DO jk = 1, jpkm1
531         DO jj = 1, jpjm1
532            DO ji = 1, fs_jpim1   ! vector opt.
533               ! i-gradient of T and S at jj
534               zwx (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk)
535               zwy (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk)
536               ! j-gradient of T and S at jj
537               zwz (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk)
538               zwt (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk)
539            END DO
540         END DO
541      END DO
542#  if defined key_partial_steps
543      ! partial steps correction at the bottom ocean level
544      DO jj = 1, jpjm1
545         DO ji = 1, fs_jpim1   ! vector opt.
546            ! last ocean level
547            iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
548            ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
549            ! i-gradient of T and S
550            zwx (ji,jj,iku) = gtu(ji,jj)
551            zwy (ji,jj,iku) = gsu(ji,jj)
552            ! j-gradient of T and S
553            zwz (ji,jj,ikv) = gtv(ji,jj) 
554            zwt (ji,jj,ikv) = gsv(ji,jj) 
555         END DO
556         END DO
557#endif
558
559
560      ! I.2 Vertical fluxes
561      ! -------------------
562
563      ! Surface and bottom vertical fluxes set to zero
564      ztfw(:,:, 1 ) = 0.e0
565      zsfw(:,:, 1 ) = 0.e0
566      ztfw(:,:,jpk) = 0.e0
567      zsfw(:,:,jpk) = 0.e0
568
569      ! interior (2=<jk=<jpk-1)
570      DO jk = 2, jpkm1
571         DO jj = 2, jpjm1
572            DO ji = fs_2, fs_jpim1   ! vector opt.
573               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
574
575               zmku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
576                              + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
577
578               zmkv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
579                              + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
580
581               zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
582               zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
583
584               ztfw(ji,jj,jk) = zcoef3 * (   zwx(ji  ,jj  ,jk-1) + zwx(ji-1,jj  ,jk)      &
585                                           + zwx(ji-1,jj  ,jk-1) + zwx(ji  ,jj  ,jk)  )   &
586                              + zcoef4 * (   zwz(ji  ,jj  ,jk-1) + zwz(ji  ,jj-1,jk)      &
587                                           + zwz(ji  ,jj-1,jk-1) + zwz(ji  ,jj  ,jk)  )
588
589               zsfw(ji,jj,jk) = zcoef3 * (   zwy(ji  ,jj  ,jk-1) + zwy(ji-1,jj  ,jk)      &
590                                           + zwy(ji-1,jj  ,jk-1) + zwy(ji  ,jj  ,jk)  )   &
591                              + zcoef4 * (   zwt(ji  ,jj  ,jk-1) + zwt(ji  ,jj-1,jk)      &
592                                           + zwt(ji  ,jj-1,jk-1) + zwt(ji  ,jj  ,jk)  )
593            END DO
594         END DO
595      END DO
596
597#if defined key_traldf_eiv
598      !                              ! ---------------------------------------!
599      !                              ! Eddy induced vertical advective fluxes !
600      !                              ! ---------------------------------------!
601         zwx(:,:, 1 ) = 0.e0
602         zwy(:,:, 1 ) = 0.e0
603         zwx(:,:,jpk) = 0.e0
604         zwy(:,:,jpk) = 0.e0
605
606         DO jk = 2, jpkm1
607            DO jj = 2, jpjm1
608               DO ji = fs_2, fs_jpim1   ! vector opt.
609#   if defined key_traldf_c2d || defined key_traldf_c3d
610                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
611                     &  * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
612                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
613                     &  * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
614                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
615                     &  * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
616                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
617                     &  * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
618
619                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
620#   else
621                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
622                     &  * e2u(ji-1,jj) * umask(ji-1,jj,jk)
623                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
624                     &  * e2u(ji  ,jj) * umask(ji  ,jj,jk)
625                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
626                     &  * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
627                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
628                     &  * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
629
630                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk)   &
631                     &            * ( zuwk - zuwki + zvwk - zvwki )
632#   endif
633                  zwx(ji,jj,jk) = + zcoeg3 * ( tb(ji,jj,jk) + tb(ji,jj,jk-1) )
634                  zwy(ji,jj,jk) = + zcoeg3 * ( sb(ji,jj,jk) + sb(ji,jj,jk-1) )
635
636                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + zwx(ji,jj,jk)
637                  zsfw(ji,jj,jk) = zsfw(ji,jj,jk) + zwy(ji,jj,jk)
638#   if defined key_diaeiv
639                  w_eiv(ji,jj,jk) = -2. * zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
640#   endif
641               END DO
642            END DO
643         END DO
644#endif
645
646      ! I.5 Divergence of vertical fluxes added to the general tracer trend
647      ! -------------------------------------------------------------------
648
649      DO jk = 1, jpkm1
650         DO jj = 2, jpjm1
651            DO ji = fs_2, fs_jpim1   ! vector opt.
652               zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
653               ztav = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
654               zsav = (  zsfw(ji,jj,jk) - zsfw(ji,jj,jk+1)  ) * zbtr
655               ta(ji,jj,jk) = ta(ji,jj,jk) + ztav
656               sa(ji,jj,jk) = sa(ji,jj,jk) + zsav
657            END DO
658         END DO
659      END DO
660
661      ! save the trends for diagnostic
662      !  WARNING jpttddoe is used here for vertical Gent velocity trend not for damping !!!
663      IF( l_trdtra )   THEN
664#   if defined key_traldf_eiv
665         ! Compute the vertical Gent velocity trend
666         DO jk = 1, jpkm1
667            DO jj = 2, jpjm1
668               DO ji = fs_2, fs_jpim1   ! vector opt.
669                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
670                  tavg(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * zbtr
671                  savg(ji,jj,jk) = ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) * zbtr
672               END DO
673            END DO
674         END DO
675
676         CALL trd_mod(tavg, savg, jpttddoe, 'TRA', kt)
677#   endif
678         ! Recompute the divergence of vertical fluxes ztav & zsav trends
679         ! computed at step 1.5 above in making the difference between the new
680         ! trend ta()/sa() and the previous one tdta()/tdsa() and substract
681         ! the vertical Gent velocity trend tavg()/savg() (zero if not used)
682         tavg(:,:,:) = ta(:,:,:) - tdta(:,:,:) - tavg(:,:,:) 
683         savg(:,:,:) = sa(:,:,:) - tdsa(:,:,:) - savg(:,:,:) 
684      ENDIF
685
686   END SUBROUTINE tra_zdf_iso
687
688#else
689   !!----------------------------------------------------------------------
690   !!   Dummy module :             NO rotation of the lateral mixing tensor
691   !!----------------------------------------------------------------------
692CONTAINS
693   SUBROUTINE tra_zdf_iso_vopt( kt )              ! empty routine
694      WRITE(*,*) 'tra_zdf_iso_vopt: You should not have seen this print! error?', kt
695   END SUBROUTINE tra_zdf_iso_vopt
696#endif
697
698   !!==============================================================================
699END MODULE trazdf_iso_vopt
Note: See TracBrowser for help on using the repository browser.