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_lf.F90 in NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN/dynldf_iso_lf.F90 @ 14801

Last change on this file since 14801 was 14801, checked in by francesca, 3 years ago

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