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

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

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