source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 4460

Last change on this file since 4460 was 4460, checked in by trackstand2, 7 years ago

Addition of dbg SUM outputs and alteration of loop limits in dynldfiso and dynldf_bilapg

  • Property svn:keywords set to Id
File size: 56.5 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  !  97-07  (G. Madec)  Original code
7   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
8   !!             -   !  2004-08  (C. Talandier) New trends organization
9   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion
10   !!----------------------------------------------------------------------
11#if defined key_ldfslp   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_ldfslp'                      slopes of the direction of mixing
14   !!----------------------------------------------------------------------
15   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
16   !!                  of the lateral diffusion using isopycnal or horizon-
17   !!                  tal s-coordinate laplacian operator.
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE ldfdyn_oce      ! ocean dynamics lateral physics
22   USE ldftra_oce      ! ocean tracer   lateral physics
23   USE zdf_oce         ! ocean vertical physics
24   USE trdmod          ! ocean dynamics trends
25   USE trdmod_oce      ! ocean variables trends
26   USE ldfslp          ! iso-neutral slopes
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE in_out_manager  ! I/O manager
29   USE lib_mpp         ! MPP library
30   USE prtctl          ! Print control
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_ldf_iso           ! called by step.F90
36   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
37
38!FTRANS zdiu zdju zdiv zdjv zdj1u zdj1v zfuw zfvw :I :z
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
41
42   !! * Control permutation of array indices
43#  include "oce_ftrans.h90"
44#  include "dom_oce_ftrans.h90"
45#  include "ldfdyn_oce_ftrans.h90"
46#  include "ldftra_oce_ftrans.h90"
47#  include "ldfslp_ftrans.h90"
48#  include "zdf_oce_ftrans.h90"
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52#  include "ldfdyn_substitute.h90"
53#  include "vectopt_loop_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   INTEGER FUNCTION dyn_ldf_iso_alloc()
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
64      !!----------------------------------------------------------------------
65      ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
66         &      zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
67         !
68      IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
69   END FUNCTION dyn_ldf_iso_alloc
70
71
72   SUBROUTINE dyn_ldf_iso( kt )
73      !!----------------------------------------------------------------------
74      !!                     ***  ROUTINE dyn_ldf_iso  ***
75      !!                       
76      !! ** Purpose :   Compute the before trend of the rotated laplacian
77      !!      operator of lateral momentum diffusion except the diagonal
78      !!      vertical term that will be computed in dynzdf module. Add it
79      !!      to the general trend of momentum equation.
80      !!
81      !! ** Method :
82      !!        The momentum lateral diffusive trend is provided by a 2nd
83      !!      order operator rotated along neutral or geopotential surfaces
84      !!      (in s-coordinates).
85      !!      It is computed using before fields (forward in time) and isopyc-
86      !!      nal or geopotential slopes computed in routine ldfslp.
87      !!      Here, u and v components are considered as 2 independent scalar
88      !!      fields. Therefore, the property of splitting divergent and rota-
89      !!      tional part of the flow of the standard, z-coordinate laplacian
90      !!      momentum diffusion is lost.
91      !!      horizontal fluxes associated with the rotated lateral mixing:
92      !!      u-component:
93      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
94      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
95      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
96      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
97      !!      v-component:
98      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
99      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
100      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
101      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
102      !!      take the horizontal divergence of the fluxes:
103      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
104      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
105      !!      Add this trend to the general trend (ua,va):
106      !!         ua = ua + diffu
107      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
108      !!      should be modified for applications others than orca_r2 (!!bug)
109      !!
110      !! ** Action :
111      !!        Update (ua,va) arrays with the before geopotential biharmonic
112      !!      mixing trend.
113      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
114      !!      of the rotated operator in dynzdf module
115      !!----------------------------------------------------------------------
116      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
117#if ! defined key_z_first
118      ! Then these workspace arrays can be left as 2D
119      USE wrk_nemo, ONLY:   zjvt  => wrk_2d_3    ! 2D workspace
120      USE wrk_nemo, ONLY:   zivf  => wrk_2d_4    ! 2D workspace
121      USE wrk_nemo, ONLY:   ziut  => wrk_2d_1 , zjuf  => wrk_2d_2
122      USE wrk_nemo, ONLY:   zdku  => wrk_2d_5 , zdkv => wrk_2d_6
123      USE wrk_nemo, ONLY:   zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8
124#endif
125      USE timing,   ONLY:   timing_start, timing_stop
126      !
127      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
128      !
129      INTEGER  ::   ji, jj, jk   ! dummy loop indices
130      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
131      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
132      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
133      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
134      REAL(wp) ::   zcof, zrecip
135#if defined key_z_first
136      REAL(wp) ::   zdku, zdk1u, zdki1u, zdk1i1u, zdkim1u, zdk1im1u
137      REAL(wp) ::         zdkj1u, zdk1j1u, zdkjm1u, zdk1jm1u
138      REAL(wp) ::   zdkv, zdk1v, zdki1v, zdk1i1v, zdkim1v, zdk1im1v
139      REAL(wp) ::         zdkj1v, zdk1j1v, zdkjm1v, zdk1jm1v
140      REAL(wp) ::   aziut, aziuti1, azjuf, azjufjm1
141      REAL(wp) ::   azivf, azivfim1, azjvtj1, azjvt
142#endif
143      !!----------------------------------------------------------------------
144
145! ziut zjvt zjuf zivf :I :I :z
146!!$#if defined key_z_first
147!!$      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ziut, zjuf, zivf, zjvt
148!!$      ALLOCATE( ziut(jpi,jpj,jpk), zjuf(jpi,jpj,jpk),  &
149!!$                zivf(jpi,jpj,jpk), zjvt(jpi,jpj,jpk) )
150!!$#endif
151      CALL timing_start('dyn_ldf_iso')
152
153      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
154         CALL ctl_stop('dyn_ldf_iso: requested workspace arrays unavailable')   ;   RETURN
155      END IF
156
157      IF( kt == nit000 ) THEN
158         IF(lwp) WRITE(numout,*)
159         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
160         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
161         !                                      ! allocate dyn_ldf_bilap arrays
162         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
163      ENDIF
164
165      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
166      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
167         !
168#if defined key_z_first
169         DO jj = 2, jpjm1       ! set the slopes of iso-level
170            DO ji = fs_2, fs_jpim1   
171               DO jk = 1, mbkmax(ji,jj) ! jpk       
172#else
173         DO jk = 1, jpk         ! set the slopes of iso-level
174            DO jj = 2, jpjm1
175               DO ji = fs_2, fs_jpim1   ! vector opt.
176#endif
177                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
178                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
179                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
180                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
181               END DO
182            END DO
183         END DO
184         ! Lateral boundary conditions on the slopes
185         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
186         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
187 
188!!bug
189         IF( kt == nit000 ) then
190            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
191               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
192         endif
193!!end
194      ENDIF
195
196      !CALL timing_start('dyn_ldf_iso_hslab')
197
198#if defined key_z_first
199
200      ! Vertical u- and v-shears at level jk and jk+1
201      ! ---------------------------------------------
202      ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
203      !                             zdkv(jk=1)=zdkv(jk=2)
204!!$      DO jj = 1, jpj, 1
205!!$         DO ji = 1, jpi, 1
206!!$
207!!$            ! jk=1 special case
208!!$            !zdk1u(ji,jj,1) = ( ub(ji,jj,1) -ub(ji,jj,2) ) * umask(ji,jj,2)
209!!$            zdk1v(ji,jj,1) = ( vb(ji,jj,1) -vb(ji,jj,2) ) * vmask(ji,jj,2)
210!!$            !zdku(ji,jj,1) = zdk1u(ji,jj,1)
211!!$            zdkv(ji,jj,1) = zdk1v(ji,jj,1)
212!!$
213!!$            DO jk = 2, jpkm1, 1
214!!$               !zdk1u(ji,jj,jk) = ( ub(ji,jj,jk) -ub(ji,jj,jk+1) ) * umask(ji,jj,jk+1)
215!!$               zdk1v(ji,jj,jk) = ( vb(ji,jj,jk) -vb(ji,jj,jk+1) ) * vmask(ji,jj,jk+1)
216!!$
217!!$               !zdku(ji,jj,jk) = ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) * umask(ji,jj,jk)
218!!$               zdkv(ji,jj,jk) = ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) * vmask(ji,jj,jk)
219!!$            END DO !jk
220!!$         END DO
221!!$      END DO
222
223!!$      !                               -----f-----
224!!$      ! Horizontal fluxes on U             | 
225!!$      ! --------------------===        t   u   t
226!!$      !                                    | 
227!!$      ! i-flux at t-point             -----f-----
228!!$
229!!$      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
230!!$         DO jj = 2, jpjm1
231!!$            DO ji = 2, jpi
232!!$               DO jk = 1, jpkm1, 1
233!!$                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
234!!$
235!!$                  zmskt = 1._wp/MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
236!!$                     &              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )
237!!$
238!!$                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5_wp * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
239!!$   
240!!$                  ziut(ji,jj,jk) = (  zabe1 * ( ub(ji,jj,jk)    - ub(ji-1,jj,jk) )   &
241!!$                     &              + zcof1 * ( zdku  + zdk1im1u     &
242!!$                     &                         +zdk1u + zdkim1u )  ) * tmask(ji,jj,jk)
243!!$               END DO
244!!$            END DO
245!!$         END DO
246!!$      ELSE                   ! other coordinate system (zco or sco) : e3t
247!!$         DO jj = 2, jpjm1
248!!$            DO ji = 2, jpi   
249!!$               DO jk = 1, jpkm1, 1
250!!$                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
251!!$
252!!$                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
253!!$                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
254!!$
255!!$                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
256!!$
257!!$                  ziut(ji,jj,jk) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
258!!$                     &              + zcof1 * ( zdku(ji,jj,jk) + zdk1u(ji-1,jj,jk)     &
259!!$                     &                         +zdk1u(ji,jj,jk) + zdku (ji-1,jj,jk) )  ) * tmask(ji,jj,jk)
260!!$
261!!$               END DO
262!!$            END DO
263!!$         END DO
264!!$      ENDIF
265
266      ! j-flux at f-point
267! BLOCKABLE(ji,jj,jk)
268! BLOCKING SIZE (0)
269!!$      DO jj = 1, jpjm1, 1
270!!$! BLOCKING SIZE (0)
271!!$         DO ji = 1, jpim1, 1
272!!$! BLOCKING SIZE (4)
273!!$            DO jk = 1, jpkm1, 1
274!!$               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
275!!$
276!!$               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
277!!$                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
278!!$
279!!$               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
280!!$
281!!$               zjuf(ji,jj,jk) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
282!!$                  &              + zcof2 * ( zdku (ji,jj+1,jk) + zdk1u(ji,jj,jk)     &
283!!$                  &                         +zdk1u(ji,jj+1,jk) + zdku (ji,jj,jk) )  ) * fmask(ji,jj,jk)
284!!$            END DO
285!!$         END DO
286!!$      END DO
287
288      !                                |   t   |
289      ! Horizontal fluxes on V         |       |
290      ! --------------------===        f---v---f
291      !                                |       |
292      !                                |   t   |
293     
294!!$      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
295!!$         DO jj = 2, jpj
296!!$            DO ji = 1, jpim1
297!!$               DO jk = 1, jpkm1, 1
298!!$
299!!$                  ! i-flux at f-point              |   t   |
300!!$                  zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
301!!$
302!!$                  zmskf = 1._wp/MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
303!!$                                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp )
304!!$
305!!$                  zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
306!!$
307!!$                  zivf(ji,jj,jk) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        &
308!!$                                    + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk)     &
309!!$                                              +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) )  ) * fmask(ji,jj,jk)
310!!$
311!!$
312!!$                  ! j-flux at t-point
313!!$                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
314!!$
315!!$                  zmskt = 1._wp/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
316!!$                                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp )
317!!$
318!!$                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
319!!$
320!!$                  zjvt(ji,jj,jk) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
321!!$                                    + zcof2 * ( zdkv(ji,jj-1,jk) + zdk1v(ji,jj,jk)     &
322!!$                                               +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) )  ) * tmask(ji,jj,jk)
323!!$               END DO
324!!$            END DO
325!!$         END DO
326!!$      ELSE                   ! other coordinate system (zco or sco) : e3t
327!!$         DO jj = 2, jpj
328!!$            DO ji = 1, jpim1
329!!$               DO jk = 1, jpkm1
330!!$
331!!$                  ! i-flux at f-point
332!!$                  zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
333!!$
334!!$                  zmskf = 1._wp/MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
335!!$                                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp )
336!!$
337!!$                  zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
338!!$
339!!$                  zivf(ji,jj,jk) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        &
340!!$                                    + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk)     &
341!!$                                               +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) )  ) * fmask(ji,jj,jk)
342!!$                  ! j-flux at t-point
343!!$                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
344!!$
345!!$                  zmskt = 1._wp/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
346!!$                                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp )
347!!$
348!!$                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
349!!$
350!!$                  zjvt(ji,jj,jk) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
351!!$                                    + zcof2 * ( zdkv (ji,jj-1,jk) + zdk1v(ji,jj,jk)     &
352!!$                                               +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) )  ) * tmask(ji,jj,jk)
353!!$               END DO
354!!$            END DO
355!!$         END DO
356!!$
357!!$      ENDIF
358
359
360
361      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
362
363         DO jj = 2, jpjm1
364            DO ji = 2, jpim1
365!DIR$ SHORTLOOP
366               ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
367               ! If condition on jk>1. This does mean that there will be out-of-bounds *reads*
368               ! when jk is 1 but that doesn't matter.
369!DIR$ SAFE_ADDRESS
370               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1, 1
371
372                  ! Vertical u- and v-shears at level jk and jk+1
373                  ! ---------------------------------------------
374                  ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
375                  !                             zdkv(jk=1)=zdkv(jk=2)
376                  zdk1u    = ( ub(ji  ,jj  ,jk) -ub(ji  ,jj  ,jk+1) ) * umask(ji  ,jj,jk+1)
377                  zdk1i1u  = ( ub(ji+1,jj  ,jk) -ub(ji+1,jj  ,jk+1) ) * umask(ji+1,jj,jk+1)
378                  zdk1j1u  = ( ub(ji  ,jj+1,jk) -ub(ji  ,jj+1,jk+1) ) * umask(ji  ,jj+1,jk+1)
379                  zdk1im1u = ( ub(ji-1,jj  ,jk) -ub(ji-1,jj  ,jk+1) ) * umask(ji-1,jj,jk+1)
380                  zdk1jm1u = ( ub(ji  ,jj-1,jk) -ub(ji  ,jj-1,jk+1) ) * umask(ji  ,jj-1,jk+1)
381                  IF(jk > 1)THEN
382                     zdku   = ( ub(ji  ,jj  ,jk-1) - ub(ji  ,jj  ,jk) ) * umask(ji,jj,jk)
383                     zdki1u = ( ub(ji+1,jj  ,jk-1) - ub(ji+1,jj  ,jk) ) * umask(ji+1,jj,jk)
384                     zdkim1u= ( ub(ji-1,jj  ,jk-1) - ub(ji-1,jj  ,jk) ) * umask(ji-1,jj,jk)
385                     zdkj1u = ( ub(ji  ,jj+1,jk-1) - ub(ji  ,jj+1,jk) ) * umask(ji,jj+1,jk)
386                     zdkjm1u= ( ub(ji  ,jj-1,jk-1) - ub(ji  ,jj-1,jk) ) * umask(ji,jj-1,jk)
387                  ELSE
388                     zdku   = zdk1u
389                     zdki1u = zdk1i1u
390                     zdkim1u= zdk1im1u
391                     zdkj1u = zdk1j1u
392                     zdkjm1u= zdk1jm1u
393                  END IF
394
395                  zdku    = zdku + zdk1u
396                  zdki1u  = zdki1u + zdk1i1u
397                  zdkim1u = zdkim1u + zdk1im1u
398                  zdkj1u  = zdkj1u +zdk1j1u
399                  zdkjm1u = zdk1jm1u + zdkjm1u
400                 
401                  ! volume elements
402                  zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
403
404                  ! horizontal component of isopycnal momentum diffusive trends
405
406                  !                               -----f-----
407                  ! Horizontal fluxes on U             | 
408                  ! --------------------===        t   u   t
409                  !                                    | 
410                  ! i-flux at t-point             -----f-----
411
412                  ! z-coordinate - partial steps : min(e3u)
413                  aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji+1,jj)) * &
414                       ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & 
415                       (1._wp/MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   &
416                               + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1._wp )) * 0.5  * &
417                               ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * &
418                               ( zdku + zdki1u )  ) * tmask(ji+1,jj,jk)
419
420                  aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * &
421                    ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & 
422                    (1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
423                            + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * &
424                    ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * &
425                    ( zdku + zdkim1u )  ) * tmask(ji,jj,jk)
426
427                  ! j-flux at f-point - identical to non-ln_zps
428
429                  azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * &
430                    ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * &
431                    (1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
432                            + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * &
433                    ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) *  &
434                    ( zdku + zdkj1u )  ) * fmask(ji,jj,jk)
435
436                  azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * &
437                    ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * &
438                    (1./MAX(  umask(ji,jj,jk  )+umask(ji,jj-1,jk+1)   &
439                            + umask(ji,jj,jk+1)+umask(ji,jj-1,jk  ), 1. )) * 0.5  * &
440                    ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) *  &
441                    ( zdku  + zdkjm1u  )  ) * fmask(ji,jj-1,jk)
442
443                  ! Second derivative (divergence) and add to the general trend
444                  ! -----------------------------------------------------------
445
446                  zuah =( &
447!                       ziut (ji+1,jj,jk) - &
448                       aziuti1 - &
449!                       ziut (ji  ,jj,jk) + &
450                       aziut + &
451!                       zjuf (ji  ,jj,jk) - &
452                       azjuf - &
453!                       zjuf (ji,jj-1,jk)   &
454                       azjufjm1 &
455                       ) / zbu
456
457                  ! add the trends to the general trends
458                  ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
459
460               END DO
461
462!DIR$ SHORTLOOP
463! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
464! If condition on jk>1. This does mean that there will be out-of-bounds *reads*
465! when jk is 1 but that doesn't matter.
466!DIR$ SAFE_ADDRESS
467               DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
468
469                  zdk1v    = ( vb(ji  ,jj  ,jk) -vb(ji  ,jj  ,jk+1) ) * vmask(ji  ,jj,jk+1)
470                  zdk1i1v  = ( vb(ji+1,jj  ,jk) -vb(ji+1,jj  ,jk+1) ) * vmask(ji+1,jj,jk+1)
471                  zdk1im1v = ( vb(ji-1,jj  ,jk) -vb(ji-1,jj  ,jk+1) ) * vmask(ji-1,jj,jk+1)
472                  zdk1j1v  = ( vb(ji  ,jj+1,jk) -vb(ji  ,jj+1,jk+1) ) * vmask(ji  ,jj+1,jk+1)
473                  zdk1jm1v = ( vb(ji  ,jj-1,jk) -vb(ji  ,jj-1,jk+1) ) * vmask(ji  ,jj-1,jk+1)
474                  IF(jk > 1)THEN
475                     zdkv   = ( vb(ji  ,jj  ,jk-1) - vb(ji  ,jj  ,jk) ) * vmask(ji,jj,jk)
476                     zdki1v = ( vb(ji+1,jj  ,jk-1) - vb(ji+1,jj  ,jk) ) * vmask(ji+1,jj,jk)
477                     zdkim1v= ( vb(ji-1,jj  ,jk-1) - vb(ji-1,jj  ,jk) ) * vmask(ji-1,jj,jk)
478                     zdkj1v = ( vb(ji  ,jj+1,jk-1) - vb(ji  ,jj+1,jk) ) * vmask(ji,jj+1,jk)
479                     zdkjm1v= ( vb(ji  ,jj-1,jk-1) - vb(ji  ,jj-1,jk) ) * vmask(ji,jj-1,jk)
480                  ELSE
481                     zdkv   = zdk1v
482                     zdki1v = zdk1i1v
483                     zdkim1v= zdk1im1v
484                     zdkj1v = zdk1j1v
485                     zdkjm1v= zdk1jm1v
486                  END IF
487                 
488                  zdkv = zdkv + zdk1v
489                  zdki1v = zdk1i1v + zdki1v
490                  zdkim1v = zdkim1v +zdk1im1v
491                  zdkj1v = zdk1j1v + zdkj1v
492                  zdkjm1v = zdkjm1v +zdk1jm1v
493
494                  !                                |   t   |
495                  ! Horizontal fluxes on V         |       |
496                  ! --------------------===        f---v---f
497                  !                                |       |
498                  !                                |   t   |
499
500                  ! i-flux at f-point - identical to non-ln_zps case
501                  azivf = (  (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / &
502                    e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        &
503                    + (- aht0 * e2f(ji,jj) / &
504                    MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
505                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * &
506                    0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * &
507                    ( zdkv + zdki1v )  ) * fmask(ji,jj,jk)
508
509                  azivfim1 = (  (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * &
510                    fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) &
511                    + (- aht0 * e2f(ji-1,jj) / &
512                    MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
513                    + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1._wp ) * 0.5_wp * &
514                    ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * &
515                    ( zdkim1v + zdkv )  ) * fmask(ji-1,jj,jk)
516
517                  ! j-flux at t-point - min(e3u) instead of e3t
518                  azjvtj1 = (  &
519                    ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * MIN( fse3v(ji,jj+1,jk), fse3v(ji,jj,jk) ) / &
520                    e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) )  &
521                    + (- aht0 * e1t(ji,jj+1) / &
522                    MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
523                    + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1._wp ) * &
524                    0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
525                    ( zdkv + zdkj1v )  &
526                         ) * tmask(ji,jj+1,jk)
527
528                  azjvt = (  &
529                    ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / &
530                    e2t(ji,jj)) * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
531                    + (- aht0 * e1t(ji,jj) /MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
532                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * 0.5_wp * &
533                    ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
534                    ( zdkjm1v + zdkv )  ) * tmask(ji,jj,jk)
535
536                  ! volume elements
537                  zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
538
539
540                  ! Second derivative (divergence) and add to the general trend
541                  ! -----------------------------------------------------------
542                  zvah =( &
543!                       zivf (ji  ,jj  ,jk) - &
544                       azivf - &
545!                       zivf (ji-1,jj  ,jk) + &
546                       azivfim1 + &
547!                       zjvt (ji  ,jj+1,jk) - &
548                       azjvtj1 - &
549!                       zjvt (ji  ,jj  ,jk)   &
550                       azjvt &
551                       ) / zbv
552
553                  ! add the trend to the general trend
554                  va (ji,jj,jk) = va (ji,jj,jk) + zvah
555               END DO
556            END DO
557         END DO
558
559      ELSE                   ! other coordinate system (zco or sco) : e3t
560
561         DO jj = 2, jpjm1
562            DO ji = 2, jpim1
563!DIR$ SHORTLOOP
564               ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
565               ! If condition on jk>1. This does mean that there will be out-of-bounds *reads*
566               ! when jk is 1 but that doesn't matter.
567!DIR$ SAFE_ADDRESS
568               DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
569
570                  ! Vertical u- and v-shears at level jk and jk+1
571                  ! ---------------------------------------------
572                  ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
573                  !                             zdkv(jk=1)=zdkv(jk=2)
574                  zdk1u    = ( ub(ji  ,jj  ,jk) -ub(ji  ,jj  ,jk+1) ) * umask(ji  ,jj,jk+1)
575                  zdk1i1u  = ( ub(ji+1,jj  ,jk) -ub(ji+1,jj  ,jk+1) ) * umask(ji+1,jj,jk+1)
576                  zdk1j1u  = ( ub(ji  ,jj+1,jk) -ub(ji  ,jj+1,jk+1) ) * umask(ji  ,jj+1,jk+1)
577                  zdk1im1u = ( ub(ji-1,jj  ,jk) -ub(ji-1,jj  ,jk+1) ) * umask(ji-1,jj,jk+1)
578                  zdk1jm1u = ( ub(ji  ,jj-1,jk) -ub(ji  ,jj-1,jk+1) ) * umask(ji  ,jj-1,jk+1)
579                  IF(jk > 1)THEN
580                     zdku   = ( ub(ji  ,jj  ,jk-1) - ub(ji  ,jj  ,jk) ) * umask(ji,jj,jk)
581                     zdki1u = ( ub(ji+1,jj  ,jk-1) - ub(ji+1,jj  ,jk) ) * umask(ji+1,jj,jk)
582                     zdkim1u= ( ub(ji-1,jj  ,jk-1) - ub(ji-1,jj  ,jk) ) * umask(ji-1,jj,jk)
583                     zdkj1u = ( ub(ji  ,jj+1,jk-1) - ub(ji  ,jj+1,jk) ) * umask(ji,jj+1,jk)
584                     zdkjm1u= ( ub(ji  ,jj-1,jk-1) - ub(ji  ,jj-1,jk) ) * umask(ji,jj-1,jk)
585                  ELSE
586                     zdku   = zdk1u
587                     zdki1u = zdk1i1u
588                     zdkim1u= zdk1im1u
589                     zdkj1u = zdk1j1u
590                     zdkjm1u= zdk1jm1u
591                  END IF
592
593                  zdku    = zdku + zdk1u
594                  zdki1u  = zdki1u + zdk1i1u
595                  zdkim1u = zdkim1u + zdk1im1u
596                  zdkj1u  = zdkj1u +zdk1j1u
597                  zdkjm1u = zdk1jm1u + zdkjm1u
598                 
599                  ! volume element
600                  zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
601
602                  ! horizontal component of isopycnal momentum diffusive trends
603
604                  !                               -----f-----
605                  ! Horizontal fluxes on U             | 
606                  ! --------------------===        t   u   t
607                  !                                    | 
608                  ! i-flux at t-point             -----f-----
609
610                  ! other coordinate system (zco or sco) : e3t
611                  aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * &
612                       ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * & 
613                       (1./MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   &
614                               + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1. )) * 0.5  * &
615                               ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * &
616                               ( zdku + zdki1u )  ) * tmask(ji+1,jj,jk)
617
618                  aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * &
619                    ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * & 
620                    (1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
621                            + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * &
622                    ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * &
623                    ( zdku + zdkim1u )  ) * tmask(ji,jj,jk)
624
625                  azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * &
626                    ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * &
627                    (1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
628                            + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * &
629                    ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) *  &
630                    ( zdku + zdkj1u )  ) * fmask(ji,jj,jk)
631
632                  azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * &
633                    ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * &
634                    (1./MAX(  umask(ji,jj,jk  )+umask(ji,jj-1,jk+1)   &
635                            + umask(ji,jj,jk+1)+umask(ji,jj-1,jk  ), 1. )) * 0.5  * &
636                    ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) *  &
637                    ( zdku  + zdkjm1u  )  ) * fmask(ji,jj-1,jk)
638
639
640                  ! Second derivative (divergence) and add to the general trend
641                  ! -----------------------------------------------------------
642                  zuah =( &
643!                       ziut (ji+1,jj,jk) - &
644                       aziuti1 - &
645!                       ziut (ji  ,jj,jk) + &
646                       aziut + &
647!                       zjuf (ji  ,jj,jk) - &
648                       azjuf - &
649!                       zjuf (ji,jj-1,jk)   &
650                       azjufjm1 &
651                       ) / zbu
652
653                  ! add the trend to the general trend
654                  ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
655
656               END DO
657
658!DIR$ SHORTLOOP
659! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
660! If condition on jk>1. This does mean that there will be out-of-bounds *reads*
661! when jk is 1 but that doesn't matter.
662!DIR$ SAFE_ADDRESS
663               DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
664
665                  zdk1v    = ( vb(ji  ,jj  ,jk) -vb(ji  ,jj  ,jk+1) ) * vmask(ji  ,jj,jk+1)
666                  zdk1i1v  = ( vb(ji+1,jj  ,jk) -vb(ji+1,jj  ,jk+1) ) * vmask(ji+1,jj,jk+1)
667                  zdk1im1v = ( vb(ji-1,jj  ,jk) -vb(ji-1,jj  ,jk+1) ) * vmask(ji-1,jj,jk+1)
668                  zdk1j1v  = ( vb(ji  ,jj+1,jk) -vb(ji  ,jj+1,jk+1) ) * vmask(ji  ,jj+1,jk+1)
669                  zdk1jm1v = ( vb(ji  ,jj-1,jk) -vb(ji  ,jj-1,jk+1) ) * vmask(ji  ,jj-1,jk+1)
670                  IF(jk > 1)THEN
671                     zdkv   = ( vb(ji  ,jj  ,jk-1) - vb(ji  ,jj  ,jk) ) * vmask(ji,jj,jk)
672                     zdki1v = ( vb(ji+1,jj  ,jk-1) - vb(ji+1,jj  ,jk) ) * vmask(ji+1,jj,jk)
673                     zdkim1v= ( vb(ji-1,jj  ,jk-1) - vb(ji-1,jj  ,jk) ) * vmask(ji-1,jj,jk)
674                     zdkj1v = ( vb(ji  ,jj+1,jk-1) - vb(ji  ,jj+1,jk) ) * vmask(ji,jj+1,jk)
675                     zdkjm1v= ( vb(ji  ,jj-1,jk-1) - vb(ji  ,jj-1,jk) ) * vmask(ji,jj-1,jk)
676                  ELSE
677                     zdkv   = zdk1v
678                     zdki1v = zdk1i1v
679                     zdkim1v= zdk1im1v
680                     zdkj1v = zdk1j1v
681                     zdkjm1v= zdk1jm1v
682                  END IF
683                 
684                  zdkv = zdkv + zdk1v
685                  zdki1v = zdk1i1v + zdki1v
686                  zdkim1v = zdkim1v +zdk1im1v
687                  zdkj1v = zdk1j1v + zdkj1v
688                  zdkjm1v = zdkjm1v +zdk1jm1v
689
690                  azivf = (  (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / &
691                    e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        &
692                    + (- aht0 * e2f(ji,jj) / &
693                    MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
694                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * &
695                    0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * &
696                    ( zdkv + zdki1v )  ) * fmask(ji,jj,jk)
697
698                  azivfim1 = (  (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * &
699                    fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) &
700                    + (- aht0 * e2f(ji-1,jj) / &
701                    MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
702                    + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1._wp ) * 0.5_wp * &
703                    ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * &
704                    ( zdkim1v + zdkv )  ) * fmask(ji-1,jj,jk)
705
706                  azjvtj1 = (  &
707                    ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / &
708                    e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) )  &
709                    + (- aht0 * e1t(ji,jj+1) / &
710                    MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
711                    + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1._wp ) * &
712                    0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
713                    ( zdkv + zdkj1v )  &
714                         ) * tmask(ji,jj+1,jk)
715
716                  azjvt = (  ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * &
717                    ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
718                    + (- aht0 * e1t(ji,jj) /MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
719                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * 0.5_wp * &
720                    ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
721                    ( zdkjm1v + zdkv )  ) * tmask(ji,jj,jk)
722
723                  ! volume elements
724                  zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
725
726                  ! Second derivative (divergence) and add to the general trend
727                  ! -----------------------------------------------------------
728                  zvah =( &
729!                       zivf (ji  ,jj  ,jk) - &
730                       azivf - &
731!                       zivf (ji-1,jj  ,jk) + &
732                       azivfim1 + &
733!                       zjvt (ji  ,jj+1,jk) - &
734                       azjvtj1 - &
735!                       zjvt (ji  ,jj  ,jk)   &
736                       azjvt &
737                       ) / zbv
738
739                  ! add the trend to the general trend
740                  va (ji,jj,jk) = va (ji,jj,jk) + zvah
741               END DO
742            END DO
743         END DO
744
745      END IF ! ln_zps
746#else
747      !                                                ! ===============
748      DO jk = 1, jpkm1                                 ! Horizontal slab
749         !                                             ! ===============
750
751         ! Vertical u- and v-shears at level jk and jk+1
752         ! ---------------------------------------------
753         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
754         !                             zdkv(jk=1)=zdkv(jk=2)
755
756         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
757         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
758
759         IF( jk == 1 ) THEN
760            zdku(:,:) = zdk1u(:,:)
761            zdkv(:,:) = zdk1v(:,:)
762         ELSE
763            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
764            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
765         ENDIF
766
767         !                               -----f-----
768         ! Horizontal fluxes on U             | 
769         ! --------------------===        t   u   t
770         !                                    | 
771         ! i-flux at t-point             -----f-----
772
773         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
774            DO jj = 2, jpjm1
775               DO ji = fs_2, jpi   ! vector opt.
776                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
777
778                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
779                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
780
781                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
782   
783                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
784                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
785                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
786               END DO
787            END DO
788         ELSE                   ! other coordinate system (zco or sco) : e3t
789            DO jj = 2, jpjm1
790               DO ji = fs_2, jpi   ! vector opt.
791                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
792
793                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
794                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
795
796                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
797
798                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
799                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
800                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
801               END DO
802            END DO
803         ENDIF
804
805         ! j-flux at f-point
806         DO jj = 1, jpjm1
807            DO ji = 1, fs_jpim1   ! vector opt.
808               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
809
810               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
811                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
812
813               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
814
815               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
816                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
817                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
818            END DO
819         END DO
820
821         !                                |   t   |
822         ! Horizontal fluxes on V         |       |
823         ! --------------------===        f---v---f
824         !                                |       |
825         ! i-flux at f-point              |   t   |
826
827         DO jj = 2, jpjm1
828            DO ji = 1, fs_jpim1   ! vector opt.
829               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
830
831               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
832                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
833
834               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
835
836               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
837                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
838                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
839            END DO
840         END DO
841
842         ! j-flux at t-point
843         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
844            DO jj = 2, jpj
845               DO ji = 1, fs_jpim1   ! vector opt.
846                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
847
848                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
849                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
850
851                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
852
853                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
854                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
855                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
856               END DO
857            END DO
858         ELSE                   ! other coordinate system (zco or sco) : e3t
859            DO jj = 2, jpj
860               DO ji = 1, fs_jpim1   ! vector opt.
861                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
862
863                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
864                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
865
866                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
867
868                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
869                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
870                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
871               END DO
872            END DO
873         ENDIF
874
875
876         ! Second derivative (divergence) and add to the general trend
877         ! -----------------------------------------------------------
878
879         DO jj = 2, jpjm1
880            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
881               ! volume elements
882               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
883               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
884               ! horizontal component of isopycnal momentum diffusive trends
885               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
886                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
887               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
888                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
889               ! add the trends to the general trends
890               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
891               va (ji,jj,jk) = va (ji,jj,jk) + zvah
892            END DO
893         END DO
894         !                                             ! ===============
895      END DO                                           !   End of slab
896      !                                                ! ===============
897#endif
898      !CALL timing_stop('dyn_ldf_iso_hslab','section')
899
900      ! print sum trends (used for debugging)
901      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
902         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
903
904      !CALL timing_start('dyn_ldf_iso_vslab')
905
906#if defined key_z_first
907      !                                                ! ===============
908      DO jj = 2, jpjm1                                 !  Vertical slab
909         !                                             ! ===============
910
911         ! I. vertical trends associated with the lateral mixing
912         ! =====================================================
913         !  (excluding the vertical flux proportional to dk[t]
914
915
916         ! I.1 horizontal momentum gradient
917         ! --------------------------------
918
919         DO ji = 2, jpi
920            DO jk = 1, mbkmax(ji,jj) ! jpk
921               ! i-gradient of u at jj
922               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
923               ! j-gradient of u and v at jj
924               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
925               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
926               ! j-gradient of u and v at jj+1
927               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
928               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
929            END DO
930         END DO
931
932         DO ji = 1, jpim1
933            DO jk = 1, mbkmax(ji,jj) ! jpk
934               ! i-gradient of v at jj
935               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
936            END DO
937         END DO
938
939
940         ! I.2 Vertical fluxes
941         ! -------------------
942
943         ! Surface and bottom vertical fluxes set to zero
944!!$         DO ji = 1, jpi
945!!$            zfuw(ji, 1 ) = 0.e0
946!!$            zfvw(ji, 1 ) = 0.e0
947!!$            zfuw(ji,jpk) = 0.e0
948!!$            zfvw(ji,jpk) = 0.e0
949!!$         END DO
950
951         ! interior (2=<jk=<jpk-1) on U field
952         DO ji = 2, jpim1
953            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
954               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
955
956               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
957               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
958
959               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
960                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
961               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
962                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
963
964               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
965               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
966               ! vertical flux on u field
967               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
968                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
969                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
970                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
971               ! update avmu (add isopycnal vertical coefficient to avmu)
972               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
973               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
974            END DO
975         END DO
976
977         ! interior (2=<jk=<jpk-1) on V field
978         DO ji = 2, jpim1
979            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
980               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
981
982               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
983               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
984
985               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
986                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
987               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
988                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
989
990               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
991               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
992               ! vertical flux on v field
993               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
994                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
995                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
996                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
997               ! update avmv (add isopycnal vertical coefficient to avmv)
998               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
999               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
1000            END DO
1001!         END DO
1002
1003
1004         ! I.3 Divergence of vertical fluxes added to the general tracer trend
1005         ! -------------------------------------------------------------------
1006!         DO ji = 2, jpim1
1007            zfuw(ji,1)             = 0.0_wp
1008            zfuw(ji,mbkmax(ji,jj)) = 0.0_wp
1009            zfvw(ji,1)             = 0.0_wp
1010            zfvw(ji,mbkmax(ji,jj)) = 0.0_wp
1011
1012            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
1013               ! volume elements
1014               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
1015               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
1016               ! part of the k-component of isopycnal momentum diffusive trends
1017               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
1018               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
1019               ! add the trends to the general trends
1020               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
1021               va(ji,jj,jk) = va(ji,jj,jk) + zvav
1022            END DO
1023         END DO
1024         !                                             ! ===============
1025      END DO                                           !   End of vertical slab
1026      !                                                ! ===============
1027#else
1028      !                                                ! ===============
1029      DO jj = 2, jpjm1                                 !  Vertical slab
1030         !                                             ! ===============
1031
1032 
1033         ! I. vertical trends associated with the lateral mixing
1034         ! =====================================================
1035         !  (excluding the vertical flux proportional to dk[t]
1036
1037
1038         ! I.1 horizontal momentum gradient
1039         ! --------------------------------
1040
1041         DO jk = 1, jpk
1042            DO ji = 2, jpi
1043               ! i-gradient of u at jj
1044               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
1045               ! j-gradient of u and v at jj
1046               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
1047               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
1048               ! j-gradient of u and v at jj+1
1049               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
1050               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
1051            END DO
1052         END DO
1053         DO jk = 1, jpk
1054            DO ji = 1, jpim1
1055               ! i-gradient of v at jj
1056               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
1057            END DO
1058         END DO
1059
1060
1061         ! I.2 Vertical fluxes
1062         ! -------------------
1063
1064         ! Surface and bottom vertical fluxes set to zero
1065         DO ji = 1, jpi
1066            zfuw(ji, 1 ) = 0.e0
1067            zfvw(ji, 1 ) = 0.e0
1068            zfuw(ji,jpk) = 0.e0
1069            zfvw(ji,jpk) = 0.e0
1070         END DO
1071
1072         ! interior (2=<jk=<jpk-1) on U field
1073         DO jk = 2, jpkm1
1074            DO ji = 2, jpim1
1075               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
1076
1077               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
1078               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
1079
1080               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
1081                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
1082               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
1083                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
1084
1085               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
1086               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
1087               ! vertical flux on u field
1088               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
1089                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
1090                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
1091                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
1092               ! update avmu (add isopycnal vertical coefficient to avmu)
1093               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
1094               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
1095            END DO
1096         END DO
1097
1098         ! interior (2=<jk=<jpk-1) on V field
1099         DO jk = 2, jpkm1
1100            DO ji = 2, jpim1
1101               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
1102
1103               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
1104               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
1105
1106               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
1107                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
1108               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
1109                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
1110
1111               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
1112               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
1113               ! vertical flux on v field
1114               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
1115                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
1116                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
1117                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
1118               ! update avmv (add isopycnal vertical coefficient to avmv)
1119               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
1120               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
1121            END DO
1122         END DO
1123
1124
1125         ! I.3 Divergence of vertical fluxes added to the general tracer trend
1126         ! -------------------------------------------------------------------
1127         DO jk = 1, jpkm1
1128            DO ji = 2, jpim1
1129               ! volume elements
1130               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
1131               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
1132               ! part of the k-component of isopycnal momentum diffusive trends
1133               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
1134               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
1135               ! add the trends to the general trends
1136               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
1137               va(ji,jj,jk) = va(ji,jj,jk) + zvav
1138            END DO
1139         END DO
1140         !                                             ! ===============
1141      END DO                                           !   End of slab
1142      !                                                ! ===============
1143#endif
1144      !CALL timing_stop('dyn_ldf_iso_vslab','section')
1145#if defined ARPDBGSUM
1146      WRITE (*,*) narea,': ARPDBG: end dyn_ldf_iso SUM(ua)= ',SUM(ua),' at step=',kt
1147#endif
1148
1149      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn_ldf_iso: failed to release workspace arrays')
1150      !
1151      CALL timing_stop('dyn_ldf_iso','section')
1152      !
1153   END SUBROUTINE dyn_ldf_iso
1154
1155# else
1156   !!----------------------------------------------------------------------
1157   !!   Dummy module                           NO rotation of mixing tensor
1158   !!----------------------------------------------------------------------
1159CONTAINS
1160   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
1161      INTEGER, INTENT(in) :: kt
1162      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
1163   END SUBROUTINE dyn_ldf_iso
1164#endif
1165
1166   !!======================================================================
1167END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.