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

source: tags/nemo_dev_x5/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90 @ 1792

Last change on this file since 1792 was 106, checked in by opalod, 20 years ago

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