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

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

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