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 tags/nemo_v1_02/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v1_02/NEMO/OPA_SRC/TRA/trazdf_iso.F90 @ 3452

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

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