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/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_iso.F90 @ 13228

Last change on this file since 13228 was 13228, checked in by smasson, 4 years ago

better e3: update with trunk@13227 see #2385

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