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

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

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