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

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

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