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

Last change on this file since 9090 was 9090, checked in by flavoni, 6 years ago

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