source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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