source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

  • Property svn:keywords set to Id
File size: 21.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 ldfslp          ! iso-neutral slopes
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE in_out_manager  ! I/O manager
27   USE lib_mpp         ! MPP library
28   USE prtctl          ! Print control
29   USE wrk_nemo        ! Memory Allocation
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 "domzgr_substitute.h90"
43#  include "ldfdyn_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
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( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
57         &      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 + ahmb0 ) e2t * e3t / e1t  di[ ub ]
85      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
86      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
87      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
88      !!      v-component:
89      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
90      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
91      !!         zjvt = ( ahmt + ahmb0 ) 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      !!        Update (ua,va) arrays with the before geopotential biharmonic
103      !!      mixing trend.
104      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
105      !!      of the rotated operator in dynzdf module
106      !!----------------------------------------------------------------------
107      !
108      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
109      !
110      INTEGER  ::   ji, jj, jk   ! dummy loop indices
111      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
112      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
113      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
114      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
115      !
116      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
117      !!----------------------------------------------------------------------
118      !
119      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso')
120      !
121      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
122      !
123      IF( kt == nit000 ) THEN
124         IF(lwp) WRITE(numout,*)
125         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
126         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
127         !                                      ! allocate dyn_ldf_bilap arrays
128         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
129      ENDIF
130
131      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
132      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
133         !
134         DO jk = 1, jpk         ! set the slopes of iso-level
135            DO jj = 2, jpjm1
136               DO ji = fs_2, fs_jpim1   ! vector opt.
137                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
138                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
139                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
140                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
141               END DO
142            END DO
143         END DO
144         ! Lateral boundary conditions on the slopes
145         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
146         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
147 
148!!bug
149         IF( kt == nit000 ) then
150            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
151               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
152         endif
153!!end
154      ENDIF
155
156      !                                                ! ===============
157      DO jk = 1, jpkm1                                 ! Horizontal slab
158         !                                             ! ===============
159
160         ! Vertical u- and v-shears at level jk and jk+1
161         ! ---------------------------------------------
162         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
163         !                             zdkv(jk=1)=zdkv(jk=2)
164
165         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
166         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
167
168         IF( jk == 1 ) THEN
169            zdku(:,:) = zdk1u(:,:)
170            zdkv(:,:) = zdk1v(:,:)
171         ELSE
172            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
173            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
174         ENDIF
175
176         !                               -----f-----
177         ! Horizontal fluxes on U             | 
178         ! --------------------===        t   u   t
179         !                                    | 
180         ! i-flux at t-point             -----f-----
181
182         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
183            DO jj = 2, jpjm1
184               DO ji = fs_2, jpi   ! vector opt.
185                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
186
187                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
188                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
189
190                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
191   
192                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
193                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
194                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
195               END DO
196            END DO
197         ELSE                   ! other coordinate system (zco or sco) : e3t
198            DO jj = 2, jpjm1
199               DO ji = fs_2, jpi   ! vector opt.
200                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
201
202                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
203                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
204
205                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
206
207                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
208                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
209                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
210               END DO
211            END DO
212         ENDIF
213
214         ! j-flux at f-point
215         DO jj = 1, jpjm1
216            DO ji = 1, fs_jpim1   ! vector opt.
217               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
218
219               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
220                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
221
222               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
223
224               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
225                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
226                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
227            END DO
228         END DO
229
230         !                                |   t   |
231         ! Horizontal fluxes on V         |       |
232         ! --------------------===        f---v---f
233         !                                |       |
234         ! i-flux at f-point              |   t   |
235
236         DO jj = 2, jpjm1
237            DO ji = 1, fs_jpim1   ! vector opt.
238               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
239
240               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
241                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
242
243               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
244
245               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
246                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
247                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
248            END DO
249         END DO
250
251         ! j-flux at t-point
252         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
253            DO jj = 2, jpj
254               DO ji = 1, fs_jpim1   ! vector opt.
255                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
256
257                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
258                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
259
260                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
261
262                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
263                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
264                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
265               END DO
266            END DO
267         ELSE                   ! other coordinate system (zco or sco) : e3t
268            DO jj = 2, jpj
269               DO ji = 1, fs_jpim1   ! vector opt.
270                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
271
272                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
273                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
274
275                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
276
277                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
278                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
279                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
280               END DO
281            END DO
282         ENDIF
283
284
285         ! Second derivative (divergence) and add to the general trend
286         ! -----------------------------------------------------------
287
288         DO jj = 2, jpjm1
289            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
290               ! volume elements
291               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
292               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
293               ! horizontal component of isopycnal momentum diffusive trends
294               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
295                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
296               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
297                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
298               ! add the trends to the general trends
299               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
300               va (ji,jj,jk) = va (ji,jj,jk) + zvah
301            END DO
302         END DO
303         !                                             ! ===============
304      END DO                                           !   End of slab
305      !                                                ! ===============
306
307      ! print sum trends (used for debugging)
308      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
309         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
310
311
312      !                                                ! ===============
313      DO jj = 2, jpjm1                                 !  Vertical slab
314         !                                             ! ===============
315
316 
317         ! I. vertical trends associated with the lateral mixing
318         ! =====================================================
319         !  (excluding the vertical flux proportional to dk[t]
320
321
322         ! I.1 horizontal momentum gradient
323         ! --------------------------------
324
325         DO jk = 1, jpk
326            DO ji = 2, jpi
327               ! i-gradient of u at jj
328               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
329               ! j-gradient of u and v at jj
330               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
331               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
332               ! j-gradient of u and v at jj+1
333               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
334               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
335            END DO
336         END DO
337         DO jk = 1, jpk
338            DO ji = 1, jpim1
339               ! i-gradient of v at jj
340               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
341            END DO
342         END DO
343
344
345         ! I.2 Vertical fluxes
346         ! -------------------
347
348         ! Surface and bottom vertical fluxes set to zero
349         DO ji = 1, jpi
350            zfuw(ji, 1 ) = 0.e0
351            zfvw(ji, 1 ) = 0.e0
352            zfuw(ji,jpk) = 0.e0
353            zfvw(ji,jpk) = 0.e0
354         END DO
355
356         ! interior (2=<jk=<jpk-1) on U field
357         DO jk = 2, jpkm1
358            DO ji = 2, jpim1
359               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
360
361               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
362               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
363
364               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
365                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
366               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
367                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
368
369               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
370               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
371               ! vertical flux on u field
372               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
373                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
374                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
375                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
376               ! update avmu (add isopycnal vertical coefficient to avmu)
377               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
378               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
379            END DO
380         END DO
381
382         ! interior (2=<jk=<jpk-1) on V field
383         DO jk = 2, jpkm1
384            DO ji = 2, jpim1
385               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
386
387               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
388               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
389
390               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
391                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
392               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
393                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
394
395               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
396               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
397               ! vertical flux on v field
398               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
399                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
400                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
401                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
402               ! update avmv (add isopycnal vertical coefficient to avmv)
403               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
404               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
405            END DO
406         END DO
407
408
409         ! I.3 Divergence of vertical fluxes added to the general tracer trend
410         ! -------------------------------------------------------------------
411         DO jk = 1, jpkm1
412            DO ji = 2, jpim1
413               ! volume elements
414               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
415               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
416               ! part of the k-component of isopycnal momentum diffusive trends
417               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
418               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
419               ! add the trends to the general trends
420               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
421               va(ji,jj,jk) = va(ji,jj,jk) + zvav
422            END DO
423         END DO
424         !                                             ! ===============
425      END DO                                           !   End of slab
426      !                                                ! ===============
427      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
428      !
429      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso')
430      !
431   END SUBROUTINE dyn_ldf_iso
432
433# else
434   !!----------------------------------------------------------------------
435   !!   Dummy module                           NO rotation of mixing tensor
436   !!----------------------------------------------------------------------
437CONTAINS
438   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
439      INTEGER, INTENT(in) :: kt
440      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
441   END SUBROUTINE dyn_ldf_iso
442#endif
443
444   !!======================================================================
445END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.