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 branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 8568

Last change on this file since 8568 was 8568, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.2 - _NONE option + remove zts + see associated wiki page
  • Property svn:keywords set to Id
File size: 20.4 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   USE timing          ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_ldf_iso           ! called by step.F90
36   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
37
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
39   
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 )
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[ ub ]
85      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
86      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
87      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ]
88      !!      v-component:
89      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ]
90      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ]
91      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
92      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
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 (ua,va):
97      !!         ua = ua + 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      !!       -(ua,va) 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      !
108      INTEGER  ::   ji, jj, jk   ! dummy loop indices
109      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
110      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
111      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4            !   -      -
112      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
113      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
114      !!----------------------------------------------------------------------
115      !
116      IF( ln_timing )   CALL timing_start('dyn_ldf_iso')
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 jk = 1, jpk         ! set the slopes of iso-level
131            DO jj = 2, jpjm1
132               DO ji = 2, jpim1
133                  uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
134                  vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
135                  wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
136                  wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
137               END DO
138            END DO
139         END DO
140         ! Lateral boundary conditions on the slopes
141         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
142         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
143         !
144       ENDIF
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(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
156         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
157
158         IF( jk == 1 ) THEN
159            zdku(:,:) = zdk1u(:,:)
160            zdkv(:,:) = zdk1v(:,:)
161         ELSE
162            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
163            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * 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_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * 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 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
181   
182                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
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_n(ji,jj,jk) * 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 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
196
197                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
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_n(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 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
213
214               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
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_n(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 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
234
235               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
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_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * 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 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
251
252                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
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_n(ji,jj,jk) * 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 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
266
267                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
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               ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
280                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
281               va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
282                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
283            END DO
284         END DO
285         !                                             ! ===============
286      END DO                                           !   End of slab
287      !                                                ! ===============
288
289      ! print sum trends (used for debugging)
290      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
291         &                       tab3d_2=va, 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) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
311               ! j-gradient of u and v at jj
312               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
313               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
314               ! j-gradient of u and v at jj+1
315               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
316               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
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) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
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 * rn_aht_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 rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
360               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_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 * rn_aht_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 rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
386               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_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               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
396               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
397            END DO
398         END DO
399         !                                             ! ===============
400      END DO                                           !   End of slab
401      !                                                ! ===============
402      !
403      IF( ln_timing )   CALL timing_stop('dyn_ldf_iso')
404      !
405   END SUBROUTINE dyn_ldf_iso
406
407   !!======================================================================
408END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.