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

source: trunk/NEMO/OPA_SRC/TRA/trazdf_iso.F90 @ 285

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

nemo_v1_bugfix_002 : CT : use ln_traldf_iso logical instead of l_traldf_iso one because it should not be dependant of the use or not of the partial steps
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.9 KB
Line 
1MODULE trazdf_iso
2   !!==============================================================================
3   !!                    ***  MODULE  trazdf_iso  ***
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  : update the tracer trend with the vertical part of
11   !!                  the isopycnal or geopotential s-coord. operator and
12   !!                  the vertical diffusion
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers variables
16   USE dom_oce         ! ocean space and time domain variables
17   USE ldfslp          ! Make iso-neutral slopes available
18   USE ldftra_oce      ! ocean active tracers: lateral physics
19   USE zdf_oce         ! ocean vertical physics
20   USE zdfddm          ! ocean vertical physics: double diffusion
21   USE trdmod          ! ocean active tracers trends
22   USE trdmod_oce      ! ocean variables trends
23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE zdfkpp          ! KPP parameterisation
26   USE prtctl          ! Print control
27
28   IMPLICIT NONE
29   PRIVATE
30
31   !! * Accessibility
32   PUBLIC tra_zdf_iso    ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "ldftra_substitute.h90"
37#  include "ldfeiv_substitute.h90"
38#  include "zdfddm_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!----------------------------------------------------------------------
41   !!  OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Header$
43   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE tra_zdf_iso( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_zdf_iso  ***
50      !!
51      !! ** Purpose :
52      !!     Compute the trend due to the vertical tracer diffusion inclu-
53      !!     ding the vertical component of lateral mixing (only for second
54      !!     order operator, for fourth order it is already computed and
55      !!     add to the general trend in traldf.F) and add it to the general
56      !!     trend of the tracer equations.
57      !!
58      !! ** Method :
59      !!         The vertical component of the lateral diffusive trends is
60      !!      provided by a 2nd order operator rotated along neural or geopo-
61      !!      tential surfaces to which an eddy induced advection can be added
62      !!      It is computed using before fields (forward in time) and isopyc-
63      !!      nal or geopotential slopes computed in routine ldfslp.
64      !!
65      !!      First part: vertical trends associated with the lateral mixing
66      !!      ==========  (excluding the vertical flux proportional to dk[t] )
67      !!      vertical fluxes associated with the rotated lateral mixing:
68      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
69      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
70      !!      save avt coef. resulting from vertical physics alone in zavt:
71      !!         zavt = avt
72      !!      update and save in zavt the vertical eddy viscosity coefficient:
73      !!         avt = avt + wslpi^2+wslj^2
74      !!      add vertical Eddy Induced advective fluxes ('lk_traldf_eiv=T):
75      !!         zftw = zftw + { di[aht e2u mi(wslpi)]
76      !!                    +dj[aht e1v mj(wslpj)] } mk(tb)
77      !!      take the horizontal divergence of the fluxes:
78      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
79      !!      Add this trend to the general trend (ta,sa):
80      !!         ta = ta + difft
81      !!
82      !!      Second part: vertical trend associated with the vertical physics
83      !!      ===========  (including the vertical flux proportional to dk[t]
84      !!                  associated with the lateral mixing, through the
85      !!                  update of avt)
86      !!      The vertical diffusion of tracers (t & s) is given by:
87      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
88      !!      It is computed using a backward time scheme, t=ta.
89      !!      Surface and bottom boundary conditions: no diffusive flux on
90      !!      both tracers (bottom, applied through the masked field avt).
91      !!      Add this trend to the general trend ta,sa :
92      !!         ta = ta + dz( avt dz(t) )
93      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
94      !!
95      !!      Third part: recover avt resulting from the vertical physics
96      !!      ==========  alone, for further diagnostics (for example to
97      !!                  compute the turbocline depth in zdfmxl.F90).
98      !!         avt = zavt
99      !!         (avs = zavs if lk_zdfddm=T )
100      !!
101      !!      'key_trdtra' defined: trend saved for futher diagnostics.
102      !!
103      !!      macro-tasked on vertical slab (jj-loop)
104      !!
105      !! ** Action :
106      !!         Update (ta,sa) arrays with the before vertical diffusion trend
107      !!         Save in (ztdta,ztdsa) arrays the trends if 'key_trdtra' defined
108      !!
109      !! History :
110      !!   7.0  !  91-11  (G. Madec)  Original code
111      !!        !  92-06  (M. Imbard)  correction on tracer trend loops
112      !!        !  96-01  (G. Madec)  statement function for e3
113      !!        !  97-05  (G. Madec)  vertical component of isopycnal
114      !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
115      !!        !  00-08  (G. Madec)  double diffusive mixing
116      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
117      !!   9.0  !  04-08  (C. Talandier) New trends organization
118      !!   9.0  !  05-06  (C. Ethe )  non-local flux in KPP vertical mixing scheme
119      !!---------------------------------------------------------------------
120      !! * Modules used
121      USE oce                ,   &
122# if defined key_zdfddm
123         zavs => va,  &  ! use va as workspace
124# endif
125         zavt => ua      ! use ua as workspace
126
127      !! * Arguments
128      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
129
130      !! * Local save
131      REAL(wp), DIMENSION(jpk), SAVE ::  &
132         z2dt
133
134      !! * Local declarations
135      INTEGER ::   ji, jj, jk                 ! dummy loop indices
136      INTEGER ::   ikst, ikenm2, ikstp1       ! temporary integers
137#if defined key_partial_steps
138      INTEGER ::   iku, ikv, ikv1             ! temporary integers
139#endif
140      REAL(wp) ::   zta, zsa
141      REAL(wp) ::   &
142         zcoef0, zcoef3,        &  ! ???
143         zcoef4, zavi,          &  ! ???
144         zbtr, zmku, zmkv,      &  !
145         ztav, zsav
146      REAL(wp), DIMENSION(jpi,jpk) ::   &
147         zwd, zws, zwi,         &  ! ???
148         zwx, zwy, zwz, zwt        ! ???
149      REAL(wp), DIMENSION(jpi,jpk) ::   &
150         ztfw, zdit, zdjt, zdj1t,   &
151         zsfw, zdis, zdjs, zdj1s
152      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
153         ztavg, zsavg,                      & ! workspace arrays
154         ztdta, ztdsa                         ! workspace arrays
155#if defined key_traldf_eiv   ||   defined key_esopa
156      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
157         ztfwg, zsfwg 
158      REAL(wp) ::                       &
159         zcoeg3,                        &
160         zuwk, zvwk,                    &
161         zuwki, zvwki
162#endif
163      !!---------------------------------------------------------------------
164      !!  OPA 8.5, LODYC-IPSL (2002)
165      !!---------------------------------------------------------------------
166
167      IF( kt == nit000 ) THEN
168         IF(lwp) WRITE(numout,*)
169         IF(lwp) WRITE(numout,*) 'tra_zdf_iso : vertical mixing (including isopycnal component)'
170         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
171#if defined key_diaeiv
172         w_eiv(:,:,:) = 0.e0
173#endif
174      ENDIF
175
176      ! 0. Local constant initialization
177      ! --------------------------------
178      ztavg(:,:,:) = 0.e0
179      zsavg(:,:,:) = 0.e0
180
181      ! time step = 2 rdttra ex
182      IF( neuler == 0 .AND. kt == nit000 ) THEN
183         z2dt(:) =  rdttra(:)              ! restarting with Euler time stepping
184      ELSEIF( kt <= nit000 + 1) THEN
185         z2dt(:) = 2. * rdttra(:)          ! leapfrog
186      ENDIF
187
188      ! Save ta and sa trends
189      IF( l_trdtra )   THEN
190         ztdta(:,:,:) = ta(:,:,:) 
191         ztdsa(:,:,:) = sa(:,:,:) 
192      ENDIF
193
194      !                                                ! ===============
195      DO jj = 2, jpjm1                                 !  Vertical slab
196         !                                             ! ===============
197
198         ! I. vertical trends associated with the lateral mixing
199         ! =====================================================
200         !  (excluding the vertical flux proportional to dk[t]
201
202
203         ! I.1 horizontal tracer gradient
204         ! ------------------------------
205
206         DO jk = 1, jpkm1
207            DO ji = 1, jpim1
208               ! i-gradient of T and S at jj
209               zdit (ji,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk)
210               zdis (ji,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk)
211               ! j-gradient of T and S at jj
212               zdjt (ji,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk)
213               zdjs (ji,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk)
214               ! j-gradient of T and S at jj+1
215               zdj1t(ji,jk) = ( tb(ji,jj,jk)-tb(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
216               zdj1s(ji,jk) = ( sb(ji,jj,jk)-sb(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
217            END DO
218         END DO
219#  if defined key_partial_steps
220         ! partial steps correction at the bottom ocean level
221         DO ji = 1, jpim1
222            ! last ocean level
223            iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
224            ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
225            ikv1 = MIN( mbathy(ji,jj), mbathy(ji  ,jj-1) ) - 1
226            ! i-gradient of T and S at jj
227            zdit (ji,iku) = gtu(ji,jj)
228            zdis (ji,iku) = gsu(ji,jj)
229            ! j-gradient of T and S at jj
230            zdjt (ji,ikv) = gtv(ji,jj) 
231            zdjs (ji,ikv) = gsv(ji,jj) 
232            ! j-gradient of T and S at jj+1
233            zdj1t(ji,ikv1)= gtv(ji,jj-1)
234            zdj1s(ji,ikv1)= gsv(ji,jj-1)
235         END DO
236#endif
237
238
239         ! I.2 Vertical fluxes
240         ! -------------------
241
242         ! Surface and bottom vertical fluxes set to zero
243         ztfw(:, 1 ) = 0.e0
244         zsfw(:, 1 ) = 0.e0
245         ztfw(:,jpk) = 0.e0
246         zsfw(:,jpk) = 0.e0
247#if defined key_traldf_eiv
248         ztfwg(:,:, 1 ) = 0.e0
249         zsfwg(:,:, 1 ) = 0.e0
250         ztfwg(:,:,jpk) = 0.e0
251         zsfwg(:,:,jpk) = 0.e0
252#endif
253
254         ! interior (2=<jk=<jpk-1)
255         DO jk = 2, jpkm1
256            DO ji = 2, jpim1
257               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
258
259               zmku = 1./MAX( umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)   &
260                  &          +umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1. )
261
262               zmkv = 1./MAX( vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)   &
263                  &          +vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1. )
264
265               zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
266               zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
267
268               ztfw(ji,jk) = zcoef3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
269                  &                    +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
270                  &        + zcoef4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
271                  &                    +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )
272
273               zsfw(ji,jk) = zcoef3 * ( zdis (ji  ,jk-1) + zdis (ji-1,jk)     &
274                  &                    +zdis (ji-1,jk-1) + zdis (ji  ,jk) )   &
275                  &        + zcoef4 * ( zdjs (ji  ,jk-1) + zdj1s(ji  ,jk)     &
276                  &                    +zdj1s(ji  ,jk-1) + zdjs (ji  ,jk) )
277            END DO
278         END DO
279
280
281         ! I.3  update and save of avt (and avs if double diffusive mixing)
282         ! ---------------------------
283
284         DO jk = 2, jpkm1
285            DO ji = 2, jpim1
286
287               zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
288                  &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
289
290               ! save avt in zavt to recover avt for mixed layer depth diag.
291               zavt(ji,jj,jk) = avt(ji,jj,jk)
292               ! add isopycnal vertical coeff. to avt
293               avt(ji,jj,jk) = avt(ji,jj,jk) + zavi
294               ! same procedure on avs if necessary
295#if defined key_zdfddm
296               ! save avs in zavs to recover avs in output files
297               zavs(ji,jj,jk) = fsavs(ji,jj,jk)
298               ! add isopycnal vertical coeff. to avs
299               fsavs(ji,jj,jk) = fsavs(ji,jj,jk) + zavi
300#endif
301            END DO
302         END DO
303
304#if defined key_traldf_eiv
305         !                              ! ---------------------------------------!
306         !                              ! Eddy induced vertical advective fluxes !
307         !                              ! ---------------------------------------!
308#if defined key_traldf_c2d || defined key_traldf_c3d
309            DO jk = 2, jpkm1
310               DO ji = 2, jpim1
311                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
312                     &  * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj)*umask(ji-1,jj,jk)
313                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
314                     &  * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj)*umask(ji  ,jj,jk)
315                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
316                     &  * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
317                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
318                     &  * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
319
320                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
321   
322                  ztfwg(ji,jj,jk) = + zcoeg3 * ( tb(ji,jj,jk) + tb(ji,jj,jk-1) )
323                  zsfwg(ji,jj,jk) = + zcoeg3 * ( sb(ji,jj,jk) + sb(ji,jj,jk-1) )
324   
325                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jj,jk)
326                  zsfw(ji,jk) = zsfw(ji,jk) + zsfwg(ji,jj,jk)
327# if defined key_diaeiv
328                  w_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
329# endif
330               END DO
331            END DO
332
333#else
334            DO jk = 2, jpkm1
335               DO ji = 2, jpim1
336                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
337                     &  * e2u(ji-1,jj)*umask(ji-1,jj,jk)
338                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
339                     &  * e2u(ji  ,jj)*umask(ji  ,jj,jk)
340                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
341                     &  * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
342                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
343                     &  * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
344   
345                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk)   &
346                     &            * ( zuwk - zuwki + zvwk - zvwki )
347
348                  ztfwg(ji,jj,jk) = + zcoeg3 * ( tb(ji,jj,jk) + tb(ji,jj,jk-1) )
349                  zsfwg(ji,jj,jk) = + zcoeg3 * ( sb(ji,jj,jk) + sb(ji,jj,jk-1) )
350
351                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jj,jk)
352                  zsfw(ji,jk) = zsfw(ji,jk) + zsfwg(ji,jj,jk)
353# if defined key_diaeiv
354                  w_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
355# endif
356               END DO
357            END DO
358#endif
359
360#endif
361
362         ! I.5 Divergence of vertical fluxes added to the general tracer trend
363         ! -------------------------------------------------------------------
364
365         DO jk = 1, jpkm1
366            DO ji = 2, jpim1
367               zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
368               ztav = (  ztfw(ji,jk) - ztfw(ji,jk+1)  ) * zbtr
369               zsav = (  zsfw(ji,jk) - zsfw(ji,jk+1)  ) * zbtr
370               ta(ji,jj,jk) = ta(ji,jj,jk) + ztav
371               sa(ji,jj,jk) = sa(ji,jj,jk) + zsav
372            END DO
373         END DO
374         !                                             ! ===============
375      END DO                                           !   End of slab
376      !                                                ! ===============
377
378      ! save the trends for diagnostic
379      !  WARNING jpttddoe is used here for vertical Gent velocity trend not for damping !!!
380      IF( l_trdtra )   THEN
381#   if defined key_traldf_eiv
382         ! Compute the vertical Gent velocity trend
383         !                                                ! ===============
384         DO jj = 2, jpjm1                                 !  Vertical slab
385            !                                             ! ===============
386            DO jk = 1, jpkm1
387               DO ji = 2, jpim1
388                  zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
389                  ztavg(ji,jj,jk) = ( ztfwg(ji,jj,jk) - ztfwg(ji,jj,jk+1) ) * zbtr
390                  zsavg(ji,jj,jk) = ( zsfwg(ji,jj,jk) - zsfwg(ji,jj,jk+1) ) * zbtr
391               END DO
392            END DO
393            !                                             ! ===============
394         END DO                                           !   End of slab
395         !                                                ! ===============
396
397         CALL trd_mod(ztavg, zsavg, jpttddoe, 'TRA', kt)
398#   endif
399         ! Recompute the divergence of vertical fluxes ztav & zsav trends
400         ! computed at step 1.5 above in making the difference between the new
401         ! trend ta()/sa() and the previous one ztdta()/ztdsa() and substract
402         ! the vertical Gent velocity trend ztavg()/zsavg() (zero if not used)
403         ztavg(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - ztavg(:,:,:) 
404         zsavg(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - zsavg(:,:,:) 
405
406         ! Save the new ta and sa trends
407         ztdta(:,:,:) = ta(:,:,:) 
408         ztdsa(:,:,:) = sa(:,:,:) 
409      ENDIF
410
411      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
412         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 1- Ta: ', mask1=tmask, &
413            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
414      ENDIF
415
416      !                                                ! ===============
417      DO jj = 2, jpjm1                                 !  Vertical slab
418         !                                             ! ===============
419
420         ! II. Vertical trend associated with the vertical physics
421         ! =======================================================
422         !     (including the vertical flux proportional to dk[t] associated
423         !      with the lateral mixing, through the avt update)
424         !     dk[ avt dk[ (t,s) ] ] diffusive trends
425
426
427         ! II.0 Matrix construction
428         ! ------------------------
429
430         ! Diagonal, inferior, superior
431         ! (including the bottom boundary condition via avt masked)
432         DO jk = 1, jpkm1
433            DO ji = 2, jpim1
434               zwi(ji,jk) = - z2dt(jk) * avt(ji,jj,jk  )   &
435                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
436               zws(ji,jk) = - z2dt(jk) * avt(ji,jj,jk+1)   &
437                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
438               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
439            END DO
440         END DO
441
442         ! Surface boudary conditions
443         DO ji = 2, jpim1
444            zwi(ji,1) = 0.e0
445            zwd(ji,1) = 1. - zws(ji,1)
446         END DO
447
448
449         ! II.1. Vertical diffusion on t
450         ! ---------------------------
451
452         ! Second member construction
453#if defined key_zdfkpp
454         ! add non-local temperature flux ( in convective case only)
455         DO jk = 1, jpkm1
456            DO ji = 2, jpim1
457               ! zrhs=right hand side
458               zwy(ji,jk) = tb(ji,jj,jk) + z2dt(jk) * ta(ji,jj,jk) &
459                  &  - z2dt(jk) * ( ghats(ji,jj,jk) * avt(ji,jj,jk) - ghats(ji,jj,jk+1) * avt(ji,jj,jk+1) ) &
460                  &               * wt0(ji,jj) / fse3t(ji,jj,jk) 
461            END DO
462         END DO
463#else
464         DO jk = 1, jpkm1
465            DO ji = 2, jpim1
466               zwy(ji,jk) = tb(ji,jj,jk) + z2dt(jk) * ta(ji,jj,jk)
467            END DO
468         END DO
469#endif
470
471         ! Matrix inversion from the first level
472         ikst = 1
473#   include "zdf.matrixsolver.h90"
474
475         ! Save the masked temperature after in ta
476         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done
477         !                  it will not be done in tranxt)
478         DO jk = 1, jpkm1
479            DO ji = 2, jpim1
480               ta(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk)
481            END DO
482         END DO
483
484
485         ! II.2 Vertical diffusion on s
486         ! ---------------------------
487
488#if defined key_zdfddm
489         ! Rebuild the Matrix as avt /= avs
490
491         ! Diagonal, inferior, superior
492         ! (including the bottom boundary condition via avs masked)
493         DO jk = 1, jpkm1
494            DO ji = 2, jpim1
495               zwi(ji,jk) = - z2dt(jk) * fsavs(ji,jj,jk  )   &
496                                       /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
497               zws(ji,jk) = - z2dt(jk) * fsavs(ji,jj,jk+1)   &
498                                       /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
499               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
500            END DO
501         END DO
502
503         ! Surface boudary conditions
504         DO ji = 2, jpim1
505            zwi(ji,1) = 0.e0
506            zwd(ji,1) = 1. - zws(ji,1)
507         END DO
508#endif
509         ! Second member construction
510
511#if defined key_zdfkpp
512         ! add non-local temperature flux ( in convective case only)
513         DO jk = 1, jpkm1
514            DO ji = 2, jpim1
515               ! zrhs=right hand side
516               zwy(ji,jk) = sb(ji,jj,jk) + z2dt(jk) * sa(ji,jj,jk) &
517                  &  - z2dt(jk) * ( ghats(ji,jj,jk) * fsavs(ji,jj,jk) - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) &
518                  &               * ws0(ji,jj) / fse3t(ji,jj,jk) 
519            END DO
520         END DO
521#else
522         DO jk = 1, jpkm1
523            DO ji = 2, jpim1
524               zwy(ji,jk) = sb(ji,jj,jk) + z2dt(jk) * sa(ji,jj,jk)
525            END DO
526         END DO
527#endif
528
529         ! Matrix inversion from the first level
530         ikst = 1
531#   include "zdf.matrixsolver.h90"
532
533         ! Save the masked salinity after in sa
534         ! (c a u t i o n:  salinity not its trend, Leap-frog scheme done
535         !                  it will not be done in tranxt)
536         DO jk = 1, jpkm1
537            DO ji = 2, jpim1
538               sa(ji,jj,jk) = zwx(ji,jk)  * tmask(ji,jj,jk)
539            END DO
540         END DO
541
542
543         ! III. recover the avt (avs) resulting from vertical physics only
544         ! ===============================================================
545
546         DO jk = 2, jpkm1
547            DO ji = 2, jpim1
548               avt(ji,jj,jk) = zavt(ji,jj,jk)
549#if defined key_zdfddm
550               fsavs(ji,jj,jk) = zavs(ji,jj,jk)
551#endif
552            END DO
553         END DO
554
555         !                                             ! ===============
556      END DO                                           !   End of slab
557      !                                                ! ===============
558
559      ! save the trends for diagnostic
560      ! compute the vertical diffusive trends in substracting the previous
561      ! trends ztdta()/ztdsa() to the new one computed via dT/dt or dS/dt
562      ! i.e. with the new temperature/salinity ta/sa computed above
563      IF( l_trdtra )   THEN
564         IF( ln_traldf_iso)   THEN
565            DO jk = 1, jpkm1
566               ztdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / z2dt(jk) ) - ztdta(:,:,jk) + ztavg(:,:,jk) 
567               ztdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / z2dt(jk) ) - ztdsa(:,:,jk) + zsavg(:,:,jk) 
568            END DO
569         ELSE
570            DO jk = 1, jpkm1
571               ztdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / z2dt(jk) ) - ztdta(:,:,jk)                             
572               ztdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / z2dt(jk) ) - ztdsa(:,:,jk)                             
573            END DO
574         ENDIF
575
576         CALL trd_mod(ztdta, ztdsa, jpttdzdf, 'TRA', kt)
577      ENDIF
578
579      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
580         CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf 2- Ta: ', mask1=tmask, &
581            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask)
582      ENDIF
583
584   END SUBROUTINE tra_zdf_iso
585
586#else
587   !!----------------------------------------------------------------------
588   !!   Dummy module               NO rotation of the lateral mixing tensor
589   !!----------------------------------------------------------------------
590CONTAINS
591   SUBROUTINE tra_zdf_iso( kt )              ! empty routine
592      WRITE(*,*) 'tra_zdf_iso: You should not have seen this print! error?', kt
593   END SUBROUTINE tra_zdf_iso
594#endif
595
596   !!==============================================================================
597END MODULE trazdf_iso
Note: See TracBrowser for help on using the repository browser.