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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynldf_iso.F90 @ 12236

Last change on this file since 12236 was 12236, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/fix_sn_cfctl_ticket2328. Fully SETTE tested
  • 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
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 "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
46   !! $Id$
47   !! Software governed by the CeCILL license (see ./LICENSE)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   INTEGER FUNCTION dyn_ldf_iso_alloc()
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
54      !!----------------------------------------------------------------------
55      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
56         &      akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
57         !
58      IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
59   END FUNCTION dyn_ldf_iso_alloc
60
61
62   SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
63      !!----------------------------------------------------------------------
64      !!                     ***  ROUTINE dyn_ldf_iso  ***
65      !!                       
66      !! ** Purpose :   Compute the before trend of the rotated laplacian
67      !!      operator of lateral momentum diffusion except the diagonal
68      !!      vertical term that will be computed in dynzdf module. Add it
69      !!      to the general trend of momentum equation.
70      !!
71      !! ** Method :
72      !!        The momentum lateral diffusive trend is provided by a 2nd
73      !!      order operator rotated along neutral or geopotential surfaces
74      !!      (in s-coordinates).
75      !!      It is computed using before fields (forward in time) and isopyc-
76      !!      nal or geopotential slopes computed in routine ldfslp.
77      !!      Here, u and v components are considered as 2 independent scalar
78      !!      fields. Therefore, the property of splitting divergent and rota-
79      !!      tional part of the flow of the standard, z-coordinate laplacian
80      !!      momentum diffusion is lost.
81      !!      horizontal fluxes associated with the rotated lateral mixing:
82      !!      u-component:
83      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ uu ]
84      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
85      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ]
86      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ]
87      !!      v-component:
88      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ]
89      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ]
90      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ]
91      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
92      !!      take the horizontal divergence of the fluxes:
93      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
94      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
95      !!      Add this trend to the general trend (uu(rhs),vv(rhs)):
96      !!         uu(rhs) = uu(rhs) + diffu
97      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
98      !!      should be modified for applications others than orca_r2 (!!bug)
99      !!
100      !! ** Action :
101      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
102      !!       -(akzu,akzv) to accompt for the diagonal vertical component
103      !!                    of the rotated operator in dynzdf module
104      !!----------------------------------------------------------------------
105      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
106      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
107      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
108      !
109      INTEGER  ::   ji, jj, jk   ! dummy loop indices
110      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
111      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
112      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      -
113      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
114      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
115      !!----------------------------------------------------------------------
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
120         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
121         !                                      ! allocate dyn_ldf_bilap arrays
122         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
123      ENDIF
124
125!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
126      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
127      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
128         !
129         DO jk = 1, jpk         ! set the slopes of iso-level
130            DO jj = 2, jpjm1
131               DO ji = 2, jpim1
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 DO
137            END DO
138         END DO
139         ! Lateral boundary conditions on the slopes
140         CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1., vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1. )
141         !
142       ENDIF
143         
144      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
145     
146      !                                                ! ===============
147      DO jk = 1, jpkm1                                 ! Horizontal slab
148         !                                             ! ===============
149
150         ! Vertical u- and v-shears at level jk and jk+1
151         ! ---------------------------------------------
152         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
153         !                             zdkv(jk=1)=zdkv(jk=2)
154
155         zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
156         zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
157
158         IF( jk == 1 ) THEN
159            zdku(:,:) = zdk1u(:,:)
160            zdkv(:,:) = zdk1v(:,:)
161         ELSE
162            zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
163            zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
164         ENDIF
165
166         !                               -----f-----
167         ! Horizontal fluxes on U             | 
168         ! --------------------===        t   u   t
169         !                                    | 
170         ! i-flux at t-point             -----f-----
171
172         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
173            DO jj = 2, jpjm1
174               DO ji = fs_2, jpi   ! vector opt.
175                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
176
177                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
178                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
179
180                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
181   
182                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
183                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
184                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
185               END DO
186            END DO
187         ELSE                   ! other coordinate system (zco or sco) : e3t
188            DO jj = 2, jpjm1
189               DO ji = fs_2, jpi   ! vector opt.
190                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
191
192                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
193                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
194
195                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
196
197                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
198                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
199                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
200               END DO
201            END DO
202         ENDIF
203
204         ! j-flux at f-point
205         DO jj = 1, jpjm1
206            DO ji = 1, fs_jpim1   ! vector opt.
207               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
208
209               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
210                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
211
212               zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
213
214               zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
215                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
216                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
217            END DO
218         END DO
219
220         !                                |   t   |
221         ! Horizontal fluxes on V         |       |
222         ! --------------------===        f---v---f
223         !                                |       |
224         ! i-flux at f-point              |   t   |
225
226         DO jj = 2, jpjm1
227            DO ji = 1, fs_jpim1   ! vector opt.
228               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
229
230               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
231                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
232
233               zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
234
235               zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
236                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
237                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
238            END DO
239         END DO
240
241         ! j-flux at t-point
242         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
243            DO jj = 2, jpj
244               DO ji = 1, fs_jpim1   ! vector opt.
245                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
246
247                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
248                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
249
250                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
251
252                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
253                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
254                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
255               END DO
256            END DO
257         ELSE                   ! other coordinate system (zco or sco) : e3t
258            DO jj = 2, jpj
259               DO ji = 1, fs_jpim1   ! vector opt.
260                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
261
262                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
263                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
264
265                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
266
267                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
268                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
269                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
270               END DO
271            END DO
272         ENDIF
273
274
275         ! Second derivative (divergence) and add to the general trend
276         ! -----------------------------------------------------------
277         DO jj = 2, jpjm1
278            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
279               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
280                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / 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) / e3v(ji,jj,jk,Kmm)
283            END DO
284         END DO
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) / e3u(ji,jj,jk,Kmm)
396               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
397            END DO
398         END DO
399         !                                             ! ===============
400      END DO                                           !   End of slab
401      !                                                ! ===============
402   END SUBROUTINE dyn_ldf_iso
403
404   !!======================================================================
405END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.