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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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