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

Last change on this file since 503 was 455, checked in by opalod, 18 years ago

nemo_v1_update_048:RB: reorganization of dynamics part, in addition change atsk to jki, suppress dynhpg_atsk.F90 dynzdf_imp_atsk.F90 dynzdf_iso.F90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
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! ajout dynzdf_iso
107      REAL(wp) ::   &
108         zcoef0, zcoef3, zcoef4, zmkt, zmkf,   &
109         zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj
110      REAL(wp), DIMENSION(jpi,jpj) ::        &
111         zfuw, zdiu, zdju, zdj1u,            & !    "        "
112         zfvw, zdiv, zdjv, zdj1v
113
114      !!----------------------------------------------------------------------
115
116      IF( kt == nit000 ) THEN
117         IF(lwp) WRITE(numout,*)
118         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
119         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
120      ENDIF
121
122!     ! s-coordinate: Iso-level diffusion on momentum but not on tracer
123      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
124 
125         ! set the slopes of iso-level
126         DO jk = 1, jpk
127            DO jj = 2, jpjm1
128               DO ji = fs_2, fs_jpim1   ! vector opt.
129                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
130                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
131                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
132                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
133               END DO
134            END DO
135         END DO
136 
137         ! Lateral boundary conditions on the slopes
138         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
139         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
140 
141!!bug
142         if( kt == nit000 ) then
143            IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
144               &                             ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj))
145         endif
146!!end
147      ENDIF
148
149      !                                                ! ===============
150      DO jk = 1, jpkm1                                 ! Horizontal slab
151         !                                             ! ===============
152
153         ! Vertical u- and v-shears at level jk and jk+1
154         ! ---------------------------------------------
155         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
156         !                             zdkv(jk=1)=zdkv(jk=2)
157
158         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
159         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
160
161         IF( jk == 1 ) THEN
162            zdku(:,:) = zdk1u(:,:)
163            zdkv(:,:) = zdk1v(:,:)
164         ELSE
165            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
166            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
167         ENDIF
168
169         !                               -----f-----
170         ! Horizontal fluxes on U             | 
171         ! --------------------===        t   u   t
172         !                                    | 
173         ! i-flux at t-point             -----f-----
174
175         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
176            DO jj = 2, jpjm1
177               DO ji = fs_2, jpi   ! vector opt.
178                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
179
180                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
181                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
182
183                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
184   
185                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
186                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
187                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
188               END DO
189            END DO
190         ELSE                   ! other coordinate system (zco or sco) : e3t
191            DO jj = 2, jpjm1
192               DO ji = fs_2, jpi   ! vector opt.
193                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
194
195                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
196                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
197
198                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
199
200                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
201                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
202                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
203               END DO
204            END DO
205         ENDIF
206
207         ! j-flux at f-point
208         DO jj = 1, jpjm1
209            DO ji = 1, fs_jpim1   ! vector opt.
210               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
211
212               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
213                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
214
215               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
216
217               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
218                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
219                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
220            END DO
221         END DO
222
223         !                                |   t   |
224         ! Horizontal fluxes on V         |       |
225         ! --------------------===        f---v---f
226         !                                |       |
227         ! i-flux at f-point              |   t   |
228
229         DO jj = 2, jpjm1
230            DO ji = 1, fs_jpim1   ! vector opt.
231               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
232
233               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
234                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
235
236               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
237
238               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
239                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
240                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
241            END DO
242         END DO
243
244         ! j-flux at t-point
245         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
246            DO jj = 2, jpj
247               DO ji = 1, fs_jpim1   ! vector opt.
248                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
249
250                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
251                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
252
253                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
254
255                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
256                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
257                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
258               END DO
259            END DO
260         ELSE                   ! other coordinate system (zco or sco) : e3t
261            DO jj = 2, jpj
262               DO ji = 1, fs_jpim1   ! vector opt.
263                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
264
265                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
266                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
267
268                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
269
270                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
271                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
272                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
273               END DO
274            END DO
275         ENDIF
276
277
278         ! Second derivative (divergence) and add to the general trend
279         ! -----------------------------------------------------------
280
281         DO jj = 2, jpjm1
282            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
283               ! volume elements
284               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
285               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
286               ! horizontal component of isopycnal momentum diffusive trends
287               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
288                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
289               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
290                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
291               ! add the trends to the general trends
292               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
293               va (ji,jj,jk) = va (ji,jj,jk) + zvah
294            END DO
295         END DO
296         !                                             ! ===============
297      END DO                                           !   End of slab
298      !                                                ! ===============
299
300      ! print sum trends (used for debugging)
301      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
302         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
303
304
305      !                                                ! ===============
306      DO jj = 2, jpjm1                                 !  Vertical slab
307         !                                             ! ===============
308
309 
310         ! I. vertical trends associated with the lateral mixing
311         ! =====================================================
312         !  (excluding the vertical flux proportional to dk[t]
313
314
315         ! I.1 horizontal momentum gradient
316         ! --------------------------------
317
318         DO jk = 1, jpk
319            DO ji = 2, jpi
320               ! i-gradient of u at jj
321               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
322               ! j-gradient of u and v at jj
323               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
324               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
325               ! j-gradient of u and v at jj+1
326               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
327               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
328            END DO
329         END DO
330         DO jk = 1, jpk
331            DO ji = 1, jpim1
332               ! i-gradient of v at jj
333               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
334            END DO
335         END DO
336
337
338         ! I.2 Vertical fluxes
339         ! -------------------
340
341         ! Surface and bottom vertical fluxes set to zero
342         DO ji = 1, jpi
343            zfuw(ji, 1 ) = 0.e0
344            zfvw(ji, 1 ) = 0.e0
345            zfuw(ji,jpk) = 0.e0
346            zfvw(ji,jpk) = 0.e0
347         END DO
348
349         ! interior (2=<jk=<jpk-1) on U field
350         DO jk = 2, jpkm1
351            DO ji = 2, jpim1
352               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
353
354               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
355               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
356
357               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
358                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
359               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
360                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
361
362               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
363               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
364               ! vertical flux on u field
365               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
366                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
367                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
368                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
369               ! update avmu (add isopycnal vertical coefficient to avmu)
370               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
371               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
372            END DO
373         END DO
374
375         ! interior (2=<jk=<jpk-1) on V field
376         DO jk = 2, jpkm1
377            DO ji = 2, jpim1
378               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
379
380               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
381               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
382
383               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
384                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
385               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
386                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
387
388               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
389               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
390               ! vertical flux on v field
391               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
392                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
393                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
394                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
395               ! update avmv (add isopycnal vertical coefficient to avmv)
396               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
397               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
398            END DO
399         END DO
400
401
402         ! I.3 Divergence of vertical fluxes added to the general tracer trend
403         ! -------------------------------------------------------------------
404         DO jk = 1, jpkm1
405            DO ji = 2, jpim1
406               ! volume elements
407               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
408               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
409               ! part of the k-component of isopycnal momentum diffusive trends
410               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
411               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
412               ! add the trends to the general trends
413               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
414               va(ji,jj,jk) = va(ji,jj,jk) + zvav
415            END DO
416         END DO
417         !                                             ! ===============
418      END DO                                           !   End of slab
419      !                                                ! ===============
420
421   END SUBROUTINE dyn_ldf_iso
422
423# else
424   !!----------------------------------------------------------------------
425   !!   Dummy module                           NO rotation of mixing tensor
426   !!----------------------------------------------------------------------
427CONTAINS
428   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
429      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
430   END SUBROUTINE dyn_ldf_iso
431#endif
432
433   !!======================================================================
434END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.