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 NEMO/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90 @ 15033

Last change on this file since 15033 was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling
  • Property svn:keywords set to Id
File size: 21.2 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator)
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 ldfslp          ! iso-neutral slopes
25   !
26   USE in_out_manager  ! I/O manager
27   USE lib_mpp         ! MPP library
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE prtctl          ! Print control
30#if defined key_loop_fusion
31   USE dynldf_iso_lf, ONLY: dyn_ldf_iso_lf   ! lateral mixing - loop fusion version (dyn_ldf_iso routine )
32#endif
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), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
41
42   !! * Substitutions
43#  include "do_loop_substitute.h90"
44#  include "domzgr_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
47   !! $Id$
48   !! Software governed by the CeCILL license (see ./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION dyn_ldf_iso_alloc()
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
55      !!----------------------------------------------------------------------
56      dyn_ldf_iso_alloc = 0
57      IF( .NOT. ALLOCATED( akzu ) ) THEN
58         ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,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      ENDIF
62   END FUNCTION dyn_ldf_iso_alloc
63
64
65   SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
66      !!----------------------------------------------------------------------
67      !!                     ***  ROUTINE dyn_ldf_iso  ***
68      !!                       
69      !! ** Purpose :   Compute the before trend of the rotated laplacian
70      !!      operator of lateral momentum diffusion except the diagonal
71      !!      vertical term that will be computed in dynzdf module. Add it
72      !!      to the general trend of momentum equation.
73      !!
74      !! ** Method :
75      !!        The momentum lateral diffusive trend is provided by a 2nd
76      !!      order operator rotated along neutral or geopotential surfaces
77      !!      (in s-coordinates).
78      !!      It is computed using before fields (forward in time) and isopyc-
79      !!      nal or geopotential slopes computed in routine ldfslp.
80      !!      Here, u and v components are considered as 2 independent scalar
81      !!      fields. Therefore, the property of splitting divergent and rota-
82      !!      tional part of the flow of the standard, z-coordinate laplacian
83      !!      momentum diffusion is lost.
84      !!      horizontal fluxes associated with the rotated lateral mixing:
85      !!      u-component:
86      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ uu ]
87      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
88      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ]
89      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ]
90      !!      v-component:
91      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ]
92      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ]
93      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ]
94      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
95      !!      take the horizontal divergence of the fluxes:
96      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
97      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
98      !!      Add this trend to the general trend (uu(rhs),vv(rhs)):
99      !!         uu(rhs) = uu(rhs) + diffu
100      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
101      !!      should be modified for applications others than orca_r2 (!!bug)
102      !!
103      !! ** Action :
104      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
105      !!       -(akzu,akzv) 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      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
110      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
111      !
112      INTEGER  ::   ji, jj, jk   ! dummy loop indices
113      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
114      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
115      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      -
116      REAL(wp), DIMENSION(A2D(nn_hls))      ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
117      REAL(wp), DIMENSION(A2D(nn_hls))      ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
118      REAL(wp), DIMENSION(A1Di(nn_hls),jpk) ::   zfuw, zdiu, zdju, zdj1u  !  -      -
119      REAL(wp), DIMENSION(A1Di(nn_hls),jpk) ::   zfvw, zdiv, zdjv, zdj1v  !  -      -
120      !!----------------------------------------------------------------------
121      !
122#if defined key_loop_fusion
123      CALL dyn_ldf_iso_lf( kt, Kbb, Kmm, puu, pvv, Krhs    )
124#else
125
126      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
127         IF( kt == nit000 ) THEN
128            IF(lwp) WRITE(numout,*)
129            IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
130            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
131            !                                      ! allocate dyn_ldf_iso arrays
132            IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
133         ENDIF
134      ENDIF
135
136!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
137      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
138      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
139         !
140         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )      ! set the slopes of iso-level
141            uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
142            vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
143            wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
144            wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
145         END_3D
146         ! Lateral boundary conditions on the slopes
147         IF (nn_hls == 1) CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
148         !
149      ENDIF
150         
151      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
152     
153      !                                                ! ===============
154      DO jk = 1, jpkm1                                 ! Horizontal slab
155         !                                             ! ===============
156
157         ! Vertical u- and v-shears at level jk and jk+1
158         ! ---------------------------------------------
159         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
160         !                             zdkv(jk=1)=zdkv(jk=2)
161
162         DO_2D( 1, 1, 1, 1 )
163            zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1)
164            zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1)
165         END_2D
166
167         IF( jk == 1 ) THEN
168            zdku(:,:) = zdk1u(:,:)
169            zdkv(:,:) = zdk1v(:,:)
170         ELSE
171            DO_2D( 1, 1, 1, 1 )
172               zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
173               zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
174            END_2D
175         ENDIF
176
177         !                               -----f-----
178         ! Horizontal fluxes on U             | 
179         ! --------------------===        t   u   t
180         !                                    | 
181         ! i-flux at t-point             -----f-----
182
183         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
184            DO_2D( 0, 1, 0, 0 )
185               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
186                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
187                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
188
189               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
190                  &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
191
192               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
193
194               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
195                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
196                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
197            END_2D
198         ELSE                   ! other coordinate system (zco or sco) : e3t
199            DO_2D( 0, 1, 0, 0 )
200               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
201                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
202
203               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
204                  &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
205
206               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
207
208               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
209                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
210                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
211            END_2D
212         ENDIF
213
214         ! j-flux at f-point
215         DO_2D( 1, 0, 1, 0 )
216            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
217               &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_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 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
223
224            zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
225               &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
226               &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
227         END_2D
228
229         !                                |   t   |
230         ! Horizontal fluxes on V         |       |
231         ! --------------------===        f---v---f
232         !                                |       |
233         ! i-flux at f-point              |   t   |
234
235         DO_2D( 1, 0, 0, 0 )
236            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
237               &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
238
239            zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
240               &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
241
242            zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
243
244            zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
245               &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
246               &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
247         END_2D
248
249         ! j-flux at t-point
250         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
251            DO_2D( 1, 0, 0, 1 )
252               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
253                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
254                  &            e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
255
256               zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
257                  &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
258
259               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
260
261               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
262                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
263                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
264            END_2D
265         ELSE                   ! other coordinate system (zco or sco) : e3t
266            DO_2D( 1, 0, 0, 1 )
267               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
268                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_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 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
274
275               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
276                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
277                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
278            END_2D
279         ENDIF
280
281
282         ! Second derivative (divergence) and add to the general trend
283         ! -----------------------------------------------------------
284         DO_2D( 0, 0, 0, 0 )      !!gm Question vectop possible??? !!bug
285            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
286               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
287               &                           / e3u(ji,jj,jk,Kmm)
288            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
289               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   &
290               &                           / e3v(ji,jj,jk,Kmm)
291         END_2D
292         !                                             ! ===============
293      END DO                                           !   End of slab
294      !                                                ! ===============
295
296      ! print sum trends (used for debugging)
297      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
298         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
299
300
301      !                                                ! ===============
302      DO jj = ntsj, ntej                               !  Vertical slab
303         !                                             ! ===============
304
305 
306         ! I. vertical trends associated with the lateral mixing
307         ! =====================================================
308         !  (excluding the vertical flux proportional to dk[t]
309
310
311         ! I.1 horizontal momentum gradient
312         ! --------------------------------
313
314         DO jk = 1, jpk
315            DO ji = ntsi, ntei + nn_hls
316               ! i-gradient of u at jj
317               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) )
318               ! j-gradient of u and v at jj
319               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) )
320               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) )
321               ! j-gradient of u and v at jj+1
322               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) )
323               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
324            END DO
325         END DO
326         DO jk = 1, jpk
327            DO ji = ntsi - nn_hls, ntei
328               ! i-gradient of v at jj
329               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
330            END DO
331         END DO
332
333
334         ! I.2 Vertical fluxes
335         ! -------------------
336
337         ! Surface and bottom vertical fluxes set to zero
338         DO ji = ntsi - nn_hls, ntei + nn_hls
339            zfuw(ji, 1 ) = 0.e0
340            zfvw(ji, 1 ) = 0.e0
341            zfuw(ji,jpk) = 0.e0
342            zfvw(ji,jpk) = 0.e0
343         END DO
344
345         ! interior (2=<jk=<jpk-1) on U field
346         DO jk = 2, jpkm1
347            DO ji = ntsi, ntei
348               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
349               !
350               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
351               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
352               !
353               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
354                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
355               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
356                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
357
358               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
359               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
360               ! vertical flux on u field
361               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
362                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
363                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
364                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
365               ! vertical mixing coefficient (akzu)
366               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
367               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
368            END DO
369         END DO
370
371         ! interior (2=<jk=<jpk-1) on V field
372         DO jk = 2, jpkm1
373            DO ji = ntsi, ntei
374               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
375               !
376               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
377               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
378               !
379               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
380                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
381               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
382                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
383
384               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
385               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
386               ! vertical flux on v field
387               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
388                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
389                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
390                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
391               ! vertical mixing coefficient (akzv)
392               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
393               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
394            END DO
395         END DO
396
397
398         ! I.3 Divergence of vertical fluxes added to the general tracer trend
399         ! -------------------------------------------------------------------
400         DO jk = 1, jpkm1
401            DO ji = ntsi, ntei
402               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   &
403                  &               / e3u(ji,jj,jk,Kmm)
404               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   &
405                  &               / e3v(ji,jj,jk,Kmm)
406            END DO
407         END DO
408         !                                             ! ===============
409      END DO                                           !   End of slab
410      !                                                ! ===============
411#endif
412   END SUBROUTINE dyn_ldf_iso
413
414   !!======================================================================
415END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.