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

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 4596

Last change on this file since 4596 was 4596, checked in by gm, 10 years ago

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