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.
dynldf_iso.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago

Add svn Id (first try), see ticket #210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 20.2 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6#if defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                      slopes of the direction of mixing
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
11   !!                  of the lateral diffusion using isopycnal or horizon-
12   !!                  tal s-coordinate laplacian operator.
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE ldfdyn_oce      ! ocean dynamics lateral physics
18   USE ldftra_oce      ! ocean tracer   lateral physics
19   USE zdf_oce         ! ocean vertical physics
20   USE trdmod          ! ocean dynamics trends
21   USE trdmod_oce      ! ocean variables trends
22   USE ldfslp          ! iso-neutral slopes
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE in_out_manager  ! I/O manager
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Routine accessibility
31   PUBLIC dyn_ldf_iso           ! called by step.F90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "ldfdyn_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   OPA 9.0 , LOCEAN-IPSL (2005)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dyn_ldf_iso( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_iso  ***
46      !!                       
47      !! ** Purpose :   Compute the before trend of the rotated laplacian
48      !!      operator of lateral momentum diffusion except the diagonal
49      !!      vertical term that will be computed in dynzdf module. Add it
50      !!      to the general trend of momentum equation.
51      !!
52      !! ** Method :
53      !!        The momentum lateral diffusive trend is provided by a 2nd
54      !!      order operator rotated along neutral or geopotential surfaces
55      !!      (in s-coordinates).
56      !!      It is computed using before fields (forward in time) and isopyc-
57      !!      nal or geopotential slopes computed in routine ldfslp.
58      !!      Here, u and v components are considered as 2 independent scalar
59      !!      fields. Therefore, the property of splitting divergent and rota-
60      !!      tional part of the flow of the standard, z-coordinate laplacian
61      !!      momentum diffusion is lost.
62      !!      horizontal fluxes associated with the rotated lateral mixing:
63      !!      u-component:
64      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
65      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
66      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
67      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
68      !!      v-component:
69      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
70      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
71      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
72      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
73      !!      take the horizontal divergence of the fluxes:
74      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
75      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
76      !!      Add this trend to the general trend (ua,va):
77      !!         ua = ua + diffu
78      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
79      !!      should be modified for applications others than orca_r2 (!!bug)
80      !!
81      !! ** Action :
82      !!        Update (ua,va) arrays with the before geopotential biharmonic
83      !!      mixing trend.
84      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
85      !!      of the rotated operator in dynzdf module
86      !!
87      !! History :
88      !!   8.0  !  97-07  (G. Madec)  Original code
89      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
90      !!   9.0  !  04-08  (C. Talandier) New trends organization
91      !!        !  05-11  (G. Madec)  s-coordinate: horizontal diffusion
92      !!----------------------------------------------------------------------
93      !! * Arguments
94      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
95
96      !! * Local declarations
97      INTEGER  ::   ji, jj, jk            ! dummy loop indices
98      REAL(wp) ::   &
99         zabe1, zabe2, zcof1, zcof2,   &  ! temporary scalars
100         zmskt, zmskf, zbu, zbv,       &
101         zuah, zvah
102      REAL(wp), DIMENSION(jpi,jpj) ::   &
103         ziut, zjuf, zjvt, zivf,        & ! temporary workspace
104         zdku, zdk1u, zdkv, zdk1v
105
106      REAL(wp) ::   &
107         zcoef0, zcoef3, zcoef4, zmkt, zmkf,   &
108         zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj
109      REAL(wp), DIMENSION(jpi,jpk) ::        &
110         zfuw, zdiu, zdju, zdj1u,            & !    "        "
111         zfvw, zdiv, zdjv, zdj1v
112
113      !!----------------------------------------------------------------------
114
115      IF( kt == nit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
119      ENDIF
120
121!     ! s-coordinate: Iso-level diffusion on momentum but not on tracer
122      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
123 
124         ! set the slopes of iso-level
125         DO jk = 1, jpk
126            DO jj = 2, jpjm1
127               DO ji = fs_2, fs_jpim1   ! vector opt.
128                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
129                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
130                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
131                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
132               END DO
133            END DO
134         END DO
135 
136         ! Lateral boundary conditions on the slopes
137         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
138         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
139 
140!!bug
141         if( kt == nit000 ) then
142            IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
143               &                             ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj))
144         endif
145!!end
146      ENDIF
147
148      !                                                ! ===============
149      DO jk = 1, jpkm1                                 ! Horizontal slab
150         !                                             ! ===============
151
152         ! Vertical u- and v-shears at level jk and jk+1
153         ! ---------------------------------------------
154         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
155         !                             zdkv(jk=1)=zdkv(jk=2)
156
157         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
158         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
159
160         IF( jk == 1 ) THEN
161            zdku(:,:) = zdk1u(:,:)
162            zdkv(:,:) = zdk1v(:,:)
163         ELSE
164            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
165            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
166         ENDIF
167
168         !                               -----f-----
169         ! Horizontal fluxes on U             | 
170         ! --------------------===        t   u   t
171         !                                    | 
172         ! i-flux at t-point             -----f-----
173
174         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
175            DO jj = 2, jpjm1
176               DO ji = fs_2, jpi   ! vector opt.
177                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
178
179                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
180                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
181
182                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
183   
184                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
185                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
186                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
187               END DO
188            END DO
189         ELSE                   ! other coordinate system (zco or sco) : e3t
190            DO jj = 2, jpjm1
191               DO ji = fs_2, jpi   ! vector opt.
192                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
193
194                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
195                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
196
197                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
198
199                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
200                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
201                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
202               END DO
203            END DO
204         ENDIF
205
206         ! j-flux at f-point
207         DO jj = 1, jpjm1
208            DO ji = 1, fs_jpim1   ! vector opt.
209               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
210
211               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
212                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
213
214               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
215
216               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
217                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
218                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
219            END DO
220         END DO
221
222         !                                |   t   |
223         ! Horizontal fluxes on V         |       |
224         ! --------------------===        f---v---f
225         !                                |       |
226         ! i-flux at f-point              |   t   |
227
228         DO jj = 2, jpjm1
229            DO ji = 1, fs_jpim1   ! vector opt.
230               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
231
232               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
233                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
234
235               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
236
237               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
238                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
239                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
240            END DO
241         END DO
242
243         ! j-flux at t-point
244         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
245            DO jj = 2, jpj
246               DO ji = 1, fs_jpim1   ! vector opt.
247                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
248
249                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
250                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
251
252                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
253
254                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
255                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
256                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
257               END DO
258            END DO
259         ELSE                   ! other coordinate system (zco or sco) : e3t
260            DO jj = 2, jpj
261               DO ji = 1, fs_jpim1   ! vector opt.
262                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
263
264                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
265                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
266
267                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
268
269                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
270                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
271                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
272               END DO
273            END DO
274         ENDIF
275
276
277         ! Second derivative (divergence) and add to the general trend
278         ! -----------------------------------------------------------
279
280         DO jj = 2, jpjm1
281            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
282               ! volume elements
283               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
284               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
285               ! horizontal component of isopycnal momentum diffusive trends
286               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
287                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
288               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
289                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
290               ! add the trends to the general trends
291               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
292               va (ji,jj,jk) = va (ji,jj,jk) + zvah
293            END DO
294         END DO
295         !                                             ! ===============
296      END DO                                           !   End of slab
297      !                                                ! ===============
298
299      ! print sum trends (used for debugging)
300      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
301         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
302
303
304      !                                                ! ===============
305      DO jj = 2, jpjm1                                 !  Vertical slab
306         !                                             ! ===============
307
308 
309         ! I. vertical trends associated with the lateral mixing
310         ! =====================================================
311         !  (excluding the vertical flux proportional to dk[t]
312
313
314         ! I.1 horizontal momentum gradient
315         ! --------------------------------
316
317         DO jk = 1, jpk
318            DO ji = 2, jpi
319               ! i-gradient of u at jj
320               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
321               ! j-gradient of u and v at jj
322               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
323               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
324               ! j-gradient of u and v at jj+1
325               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
326               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
327            END DO
328         END DO
329         DO jk = 1, jpk
330            DO ji = 1, jpim1
331               ! i-gradient of v at jj
332               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
333            END DO
334         END DO
335
336
337         ! I.2 Vertical fluxes
338         ! -------------------
339
340         ! Surface and bottom vertical fluxes set to zero
341         DO ji = 1, jpi
342            zfuw(ji, 1 ) = 0.e0
343            zfvw(ji, 1 ) = 0.e0
344            zfuw(ji,jpk) = 0.e0
345            zfvw(ji,jpk) = 0.e0
346         END DO
347
348         ! interior (2=<jk=<jpk-1) on U field
349         DO jk = 2, jpkm1
350            DO ji = 2, jpim1
351               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
352
353               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
354               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
355
356               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
357                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
358               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
359                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
360
361               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
362               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
363               ! vertical flux on u field
364               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
365                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
366                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
367                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
368               ! update avmu (add isopycnal vertical coefficient to avmu)
369               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
370               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
371            END DO
372         END DO
373
374         ! interior (2=<jk=<jpk-1) on V field
375         DO jk = 2, jpkm1
376            DO ji = 2, jpim1
377               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
378
379               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
380               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
381
382               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
383                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
384               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
385                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
386
387               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
388               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
389               ! vertical flux on v field
390               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
391                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
392                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
393                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
394               ! update avmv (add isopycnal vertical coefficient to avmv)
395               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
396               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
397            END DO
398         END DO
399
400
401         ! I.3 Divergence of vertical fluxes added to the general tracer trend
402         ! -------------------------------------------------------------------
403         DO jk = 1, jpkm1
404            DO ji = 2, jpim1
405               ! volume elements
406               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
407               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
408               ! part of the k-component of isopycnal momentum diffusive trends
409               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
410               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
411               ! add the trends to the general trends
412               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
413               va(ji,jj,jk) = va(ji,jj,jk) + zvav
414            END DO
415         END DO
416         !                                             ! ===============
417      END DO                                           !   End of slab
418      !                                                ! ===============
419
420   END SUBROUTINE dyn_ldf_iso
421
422# else
423   !!----------------------------------------------------------------------
424   !!   Dummy module                           NO rotation of mixing tensor
425   !!----------------------------------------------------------------------
426CONTAINS
427   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
428      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
429   END SUBROUTINE dyn_ldf_iso
430#endif
431
432   !!======================================================================
433END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.