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

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

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