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

source: branches/UKMO/CO6_shelfclimate_fabm/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 11738

Last change on this file since 11738 was 7567, checked in by hadjt, 7 years ago

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