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/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.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

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