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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 5760

Last change on this file since 5760 was 5758, checked in by gm, 9 years ago

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