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

Last change on this file since 34 was 34, checked in by opalod, 20 years ago

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