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

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

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