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

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

CT : BUGFIX017 : 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
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 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,        &  ! ???
128         zcoef4, zavi,          &  ! ???
129         zbtr, zmku, zmkv,      &  !
130         ztav, zsav
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,   &
137         zavt
138#if defined key_traldf_eiv   ||   defined key_esopa
139      REAL(wp), DIMENSION(jpi,jpk) ::   &
140         ztfwg, zsfwg, zavs 
141      REAL(wp) ::                       &
142         zcoeg3,                        &
143         zuwk, zvwk,                    &
144         zuwki, zvwki
145#endif
146      !!---------------------------------------------------------------------
147      !!  OPA 8.5, LODYC-IPSL (2002)
148      !!---------------------------------------------------------------------
149
150      IF( kt == nit000 ) THEN
151         IF(lwp) WRITE(numout,*)
152         IF(lwp) WRITE(numout,*) 'tra_zdf_iso : vertical mixing (including isopycnal component)'
153         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
154#if defined key_diaeiv
155         w_eiv(:,:,:) = 0.e0
156#endif
157      ENDIF
158
159      ! 0. Local constant initialization
160      ! --------------------------------
161
162      ztavg = 0.e0
163      zsavg = 0.e0
164      ! time step = 2 rdttra ex
165      IF( neuler == 0 .AND. kt == nit000 ) THEN
166         z2dt(:) =  rdttra(:)              ! restarting with Euler time stepping
167      ELSEIF( kt <= nit000 + 1) THEN
168         z2dt(:) = 2. * rdttra(:)          ! leapfrog
169      ENDIF
170
171      !                                                ! ===============
172      DO jj = 2, jpjm1                                 !  Vertical slab
173         !                                             ! ===============
174
175         ! I. vertical trends associated with the lateral mixing
176         ! =====================================================
177         !  (excluding the vertical flux proportional to dk[t]
178
179
180         ! I.1 horizontal tracer gradient
181         ! ------------------------------
182
183         DO jk = 1, jpkm1
184            DO ji = 1, jpim1
185               ! i-gradient of T and S at jj
186               zdit (ji,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk)
187               zdis (ji,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk)
188               ! j-gradient of T and S at jj
189               zdjt (ji,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk)
190               zdjs (ji,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk)
191               ! j-gradient of T and S at jj+1
192               zdj1t(ji,jk) = ( tb(ji,jj,jk)-tb(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
193               zdj1s(ji,jk) = ( sb(ji,jj,jk)-sb(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
194            END DO
195         END DO
196#  if defined key_partial_steps
197         ! partial steps correction at the bottom ocean level
198         DO ji = 1, jpim1
199            ! last ocean level
200            iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
201            ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
202            ikv1 = MIN( mbathy(ji,jj), mbathy(ji  ,jj-1) ) - 1
203            ! i-gradient of T and S at jj
204            zdit (ji,iku) = gtu(ji,jj)
205            zdis (ji,iku) = gsu(ji,jj)
206            ! j-gradient of T and S at jj
207            zdjt (ji,ikv) = gtv(ji,jj) 
208            zdjs (ji,ikv) = gsv(ji,jj) 
209            ! j-gradient of T and S at jj+1
210            zdj1t(ji,ikv1)= gtv(ji,jj-1)
211            zdj1s(ji,ikv1)= gsv(ji,jj-1)
212         END DO
213#endif
214
215
216         ! I.2 Vertical fluxes
217         ! -------------------
218
219         ! Surface and bottom vertical fluxes set to zero
220         ztfw(:, 1 ) = 0.e0
221         zsfw(:, 1 ) = 0.e0
222         ztfw(:,jpk) = 0.e0
223         zsfw(:,jpk) = 0.e0
224#if defined key_traldf_eiv
225         ztfwg(:, 1 ) = 0.e0
226         zsfwg(:, 1 ) = 0.e0
227         ztfwg(:,jpk) = 0.e0
228         zsfwg(:,jpk) = 0.e0
229#endif
230
231         ! interior (2=<jk=<jpk-1)
232         DO jk = 2, jpkm1
233            DO ji = 2, jpim1
234               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
235
236               zmku = 1./MAX( umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)   &
237                  &          +umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1. )
238
239               zmkv = 1./MAX( vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)   &
240                  &          +vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1. )
241
242               zcoef3 = zcoef0 * e2t(ji,jj) * zmku * wslpi (ji,jj,jk)
243               zcoef4 = zcoef0 * e1t(ji,jj) * zmkv * wslpj (ji,jj,jk)
244
245               ztfw(ji,jk) = zcoef3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
246                  &                    +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
247                  &        + zcoef4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
248                  &                    +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )
249
250               zsfw(ji,jk) = zcoef3 * ( zdis (ji  ,jk-1) + zdis (ji-1,jk)     &
251                  &                    +zdis (ji-1,jk-1) + zdis (ji  ,jk) )   &
252                  &        + zcoef4 * ( zdjs (ji  ,jk-1) + zdj1s(ji  ,jk)     &
253                  &                    +zdj1s(ji  ,jk-1) + zdjs (ji  ,jk) )
254            END DO
255         END DO
256
257
258         ! I.3  update and save of avt (and avs if double diffusive mixing)
259         ! ---------------------------
260
261         DO jk = 2, jpkm1
262            DO ji = 2, jpim1
263
264               zavi = fsahtw(ji,jj,jk)*( wslpi(ji,jj,jk)*wslpi(ji,jj,jk)   &
265                  &                     +wslpj(ji,jj,jk)*wslpj(ji,jj,jk) )
266
267               ! save avt in zavt to recover avt for mixed layer depth diag.
268               zavt(ji,jk) = avt(ji,jj,jk)
269               ! add isopycnal vertical coeff. to avt
270               avt(ji,jj,jk) = avt(ji,jj,jk) + zavi
271               ! same procedure on avs if necessary
272#if defined key_zdfddm
273               ! save avs in zavs to recover avs in output files
274               zavs(ji,jk) = fsavs(ji,jj,jk)
275               ! add isopycnal vertical coeff. to avs
276               fsavs(ji,jj,jk) = fsavs(ji,jj,jk) + zavi
277#endif
278            END DO
279         END DO
280
281#if defined key_traldf_eiv
282         !                              ! ---------------------------------------!
283         !                              ! Eddy induced vertical advective fluxes !
284         !                              ! ---------------------------------------!
285#if defined key_traldf_c2d || defined key_traldf_c3d
286            DO jk = 2, jpkm1
287               DO ji = 2, jpim1
288                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
289                     &  * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj)*umask(ji-1,jj,jk)
290                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
291                     &  * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj)*umask(ji  ,jj,jk)
292                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
293                     &  * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
294                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
295                     &  * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
296
297                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
298   
299                  ztfwg(ji,jk) = + zcoeg3 * ( tb(ji,jj,jk) + tb(ji,jj,jk-1) )
300                  zsfwg(ji,jk) = + zcoeg3 * ( sb(ji,jj,jk) + sb(ji,jj,jk-1) )
301   
302                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
303                  zsfw(ji,jk) = zsfw(ji,jk) + zsfwg(ji,jk)
304# if defined key_diaeiv
305                  w_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
306# endif
307               END DO
308            END DO
309
310#else
311            DO jk = 2, jpkm1
312               DO ji = 2, jpim1
313                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) )   &
314                     &  * e2u(ji-1,jj)*umask(ji-1,jj,jk)
315                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) )   &
316                     &  * e2u(ji  ,jj)*umask(ji  ,jj,jk)
317                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) )   &
318                     &  * e1v(ji,jj-1)*vmask(ji,jj-1,jk)
319                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) )   &
320                     &  * e1v(ji  ,jj)*vmask(ji  ,jj,jk)
321   
322                  zcoeg3 = + 0.25 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk)   &
323                     &            * ( zuwk - zuwki + zvwk - zvwki )
324
325                  ztfwg(ji,jk) = + zcoeg3 * ( tb(ji,jj,jk) + tb(ji,jj,jk-1) )
326                  zsfwg(ji,jk) = + zcoeg3 * ( sb(ji,jj,jk) + sb(ji,jj,jk-1) )
327
328                  ztfw(ji,jk) = ztfw(ji,jk) + ztfwg(ji,jk)
329                  zsfw(ji,jk) = zsfw(ji,jk) + zsfwg(ji,jk)
330# if defined key_diaeiv
331                  w_eiv(ji,jj,jk) = -2. *  zcoeg3 / ( e1t(ji,jj)*e2t(ji,jj) )
332# endif
333               END DO
334            END DO
335#endif
336
337#endif
338
339         ! I.5 Divergence of vertical fluxes added to the general tracer trend
340         ! -------------------------------------------------------------------
341
342         DO jk = 1, jpkm1
343            DO ji = 2, jpim1
344               zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
345               ztav = (  ztfw(ji,jk) - ztfw(ji,jk+1)  ) * zbtr
346               zsav = (  zsfw(ji,jk) - zsfw(ji,jk+1)  ) * zbtr
347               ta(ji,jj,jk) = ta(ji,jj,jk) + ztav
348               sa(ji,jj,jk) = sa(ji,jj,jk) + zsav
349#if defined key_trdtra || defined key_trdmld
350#   if defined key_traldf_eiv
351               ztavg = ( ztfwg(ji,jk) - ztfwg(ji,jk+1) ) * zbtr
352               zsavg = ( zsfwg(ji,jk) - zsfwg(ji,jk+1) ) * zbtr
353               !  WARNING ttrd(ji,jj,jk,6) used for vertical gent velocity trend
354               !                           not for damping !!!
355               ttrd(ji,jj,jk,6) = ztavg
356               strd(ji,jj,jk,6) = zsavg
357#   endif
358               ttrd(ji,jj,jk,4) = ztav - ztavg
359               strd(ji,jj,jk,4) = zsav - zsavg
360#endif
361            END DO
362         END DO
363         !                                             ! ===============
364      END DO                                           !   End of slab
365      !                                                ! ===============
366
367
368      IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
369         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
370         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
371         WRITE(numout,*) ' zdf 1- Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
372         t_ctl = zta   ;   s_ctl = zsa
373      ENDIF
374
375      !                                                ! ===============
376      DO jj = 2, jpjm1                                 !  Vertical slab
377         !                                             ! ===============
378
379         ! II. Vertical trend associated with the vertical physics
380         ! =======================================================
381         !     (including the vertical flux proportional to dk[t] associated
382         !      with the lateral mixing, through the avt update)
383         !     dk[ avt dk[ (t,s) ] ] diffusive trends
384
385
386         ! II.0 Matrix construction
387         ! ------------------------
388
389         ! Diagonal, inferior, superior
390         ! (including the bottom boundary condition via avt masked)
391         DO jk = 1, jpkm1
392            DO ji = 2, jpim1
393               zwi(ji,jk) = - z2dt(jk) * avt(ji,jj,jk  )   &
394                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
395               zws(ji,jk) = - z2dt(jk) * avt(ji,jj,jk+1)   &
396                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
397               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
398            END DO
399         END DO
400
401         ! Surface boudary conditions
402         DO ji = 2, jpim1
403            zwi(ji,1) = 0.e0
404            zwd(ji,1) = 1. - zws(ji,1)
405         END DO
406
407
408         ! II.1. Vertical diffusion on t
409         ! ---------------------------
410
411         ! Second member construction
412         DO jk = 1, jpkm1
413            DO ji = 2, jpim1
414               zwy(ji,jk) = tb(ji,jj,jk) + z2dt(jk) * ta(ji,jj,jk)
415            END DO
416         END DO
417
418         ! Matrix inversion from the first level
419         ikst = 1
420#   include "zdf.matrixsolver.h90"
421
422#if defined key_trdtra   ||   defined key_trdmld
423         ! Compute and save the vertical diffusive temperature trends
424         IF( l_traldf_iso ) THEN
425            DO jk = 1, jpkm1
426               DO ji = 2, jpim1
427                  zta = ( zwx(ji,jk) - tb(ji,jj,jk) ) / z2dt(jk)
428                  ttrd(ji,jj,jk,4) = zta - ta(ji,jj,jk) + ttrd(ji,jj,jk,4)
429               END DO
430            END DO
431         ELSE
432            DO jk = 1, jpkm1
433               DO ji = 2, jpim1
434                  zta = ( zwx(ji,jk) - tb(ji,jj,jk) ) / z2dt(jk)
435                  ttrd(ji,jj,jk,4) = zta - ta(ji,jj,jk)
436               END DO
437            END DO
438         ENDIF
439#endif
440
441         ! Save the masked temperature after in ta
442         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done
443         !                  it will not be done in tranxt)
444         DO jk = 1, jpkm1
445            DO ji = 2, jpim1
446               ta(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk)
447            END DO
448         END DO
449
450
451         ! II.2 Vertical diffusion on s
452         ! ---------------------------
453
454#if defined key_zdfddm
455         ! Rebuild the Matrix as avt /= avs
456
457         ! Diagonal, inferior, superior
458         ! (including the bottom boundary condition via avs masked)
459         DO jk = 1, jpkm1
460            DO ji = 2, jpim1
461               zwi(ji,jk) = - z2dt(jk) * fsavs(ji,jj,jk  )   &
462                                       /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
463               zws(ji,jk) = - z2dt(jk) * fsavs(ji,jj,jk+1)   &
464                                       /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
465               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
466            END DO
467         END DO
468
469         ! Surface boudary conditions
470         DO ji = 2, jpim1
471            zwi(ji,1) = 0.e0
472            zwd(ji,1) = 1. - zws(ji,1)
473         END DO
474#endif
475         ! Second member construction
476         DO jk = 1, jpkm1
477            DO ji = 2, jpim1
478               zwy(ji,jk) = sb(ji,jj,jk) + z2dt(jk) * sa(ji,jj,jk)
479            END DO
480         END DO
481
482         ! Matrix inversion from the first level
483         ikst = 1
484#   include "zdf.matrixsolver.h90"
485
486#if defined key_trdtra   ||   defined key_trdmld
487         ! Compute and save the vertical diffusive salinity trends
488         IF( l_traldf_iso ) THEN
489            DO jk = 1, jpkm1
490               DO ji = 2, jpim1
491                  zsa = ( zwx(ji,jk) - sb(ji,jj,jk) ) / z2dt(jk)
492                  strd(ji,jj,jk,4) = zsa - sa(ji,jj,jk) + strd(ji,jj,jk,4)
493               END DO
494            END DO
495         ELSE
496            DO jk = 1, jpkm1
497               DO ji = 2, jpim1
498                  zsa = ( zwx(ji,jk) - sb(ji,jj,jk) ) / z2dt(jk)
499                  strd(ji,jj,jk,4) = zsa - sa(ji,jj,jk)
500               END DO
501            END DO
502         ENDIF
503#endif
504
505         ! Save the masked salinity after in sa
506         ! (c a u t i o n:  salinity not its trend, Leap-frog scheme done
507         !                  it will not be done in tranxt)
508         DO jk = 1, jpkm1
509            DO ji = 2, jpim1
510               sa(ji,jj,jk) = zwx(ji,jk)  * tmask(ji,jj,jk)
511            END DO
512         END DO
513
514
515         ! III. recover the avt (avs) resulting from vertical physics only
516         ! ===============================================================
517
518         DO jk = 2, jpkm1
519            DO ji = 2, jpim1
520               avt(ji,jj,jk) = zavt(ji,jk)
521#if defined key_zdfddm
522               fsavs(ji,jj,jk) = zavs(ji,jk)
523#endif
524            END DO
525         END DO
526
527         !                                             ! ===============
528      END DO                                           !   End of slab
529      !                                                ! ===============
530
531      IF( l_ctl .AND. lwp ) THEN         ! print mean trends (used for debugging)
532         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
533         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
534         WRITE(numout,*) ' zdf 2- Ta: ', zta, ' Sa: ', zsa
535      ENDIF
536
537   END SUBROUTINE tra_zdf_iso
538
539#else
540   !!----------------------------------------------------------------------
541   !!   Dummy module               NO rotation of the lateral mixing tensor
542   !!----------------------------------------------------------------------
543CONTAINS
544   SUBROUTINE tra_zdf_iso( kt )              ! empty routine
545      WRITE(*,*) 'tra_zdf_iso: You should not have seen this print! error?', kt
546   END SUBROUTINE tra_zdf_iso
547#endif
548
549   !!==============================================================================
550END MODULE trazdf_iso
Note: See TracBrowser for help on using the repository browser.