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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 9441

Last change on this file since 9441 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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