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

source: tags/nemo_v1_02/NEMO/OPA_SRC/TRA/trazdf_iso_vopt.F90 @ 3452

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

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