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