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/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 10762

Last change on this file since 10762 was 10762, checked in by jcastill, 5 years ago

Revert previous changes as the removal of keywords was not uncoupled of the actual changes

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   USE wrk_nemo        ! Memory Allocation
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   dyn_ldf_iso           ! called by step.F90
37   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
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/OPA 3.3 , NEMO Consortium (2011)
46   !! $Id: dynldf_iso.F90 6140 2015-12-21 11:35:23Z timgraham $
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   INTEGER FUNCTION dyn_ldf_iso_alloc()
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
54      !!----------------------------------------------------------------------
55      ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
56         &      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 )
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[ ub ]
84      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
85      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
86      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ]
87      !!      v-component:
88      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ]
89      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ]
90      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
91      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
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 (ua,va):
96      !!         ua = ua + 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      !!        Update (ua,va) arrays with the before geopotential biharmonic
102      !!      mixing trend.
103      !!        Update (avmu,avmv) 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, zabe2, zcof1, zcof2                       ! local scalars
110      REAL(wp) ::   zmskt, zmskf                                     !   -      -
111      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
112      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
113      !
114      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
115      !!----------------------------------------------------------------------
116      !
117      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso')
118      !
119      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
120      !
121      IF( kt == nit000 ) THEN
122         IF(lwp) WRITE(numout,*)
123         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
124         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
125         !                                      ! allocate dyn_ldf_bilap arrays
126         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
127      ENDIF
128
129!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
130      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
131      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
132         !
133         DO jk = 1, jpk         ! set the slopes of iso-level
134            DO jj = 2, jpjm1
135               DO ji = 2, jpim1
136                  uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
137                  vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
138                  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
139                  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
140               END DO
141            END DO
142         END DO
143         ! Lateral boundary conditions on the slopes
144         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
145         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
146 
147!!bug
148         IF( kt == nit000 ) then
149            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
150               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
151         endif
152!!end
153      ENDIF
154
155      !                                                ! ===============
156      DO jk = 1, jpkm1                                 ! Horizontal slab
157         !                                             ! ===============
158
159         ! Vertical u- and v-shears at level jk and jk+1
160         ! ---------------------------------------------
161         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
162         !                             zdkv(jk=1)=zdkv(jk=2)
163
164         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
165         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
166
167         IF( jk == 1 ) THEN
168            zdku(:,:) = zdk1u(:,:)
169            zdkv(:,:) = zdk1v(:,:)
170         ELSE
171            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
172            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
173         ENDIF
174
175         !                               -----f-----
176         ! Horizontal fluxes on U             | 
177         ! --------------------===        t   u   t
178         !                                    | 
179         ! i-flux at t-point             -----f-----
180
181         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
182            DO jj = 2, jpjm1
183               DO ji = fs_2, jpi   ! vector opt.
184                  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)
185
186                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
187                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
188
189                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
190   
191                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
192                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
193                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
194               END DO
195            END DO
196         ELSE                   ! other coordinate system (zco or sco) : e3t
197            DO jj = 2, jpjm1
198               DO ji = fs_2, jpi   ! vector opt.
199                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj)
200
201                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
202                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
203
204                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
205
206                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
207                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
208                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
209               END DO
210            END DO
211         ENDIF
212
213         ! j-flux at f-point
214         DO jj = 1, jpjm1
215            DO ji = 1, fs_jpim1   ! vector opt.
216               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj)
217
218               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
219                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
220
221               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
222
223               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
224                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
225                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
226            END DO
227         END DO
228
229         !                                |   t   |
230         ! Horizontal fluxes on V         |       |
231         ! --------------------===        f---v---f
232         !                                |       |
233         ! i-flux at f-point              |   t   |
234
235         DO jj = 2, jpjm1
236            DO ji = 1, fs_jpim1   ! vector opt.
237               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj)
238
239               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
240                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
241
242               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
243
244               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
245                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
246                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
247            END DO
248         END DO
249
250         ! j-flux at t-point
251         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
252            DO jj = 2, jpj
253               DO ji = 1, fs_jpim1   ! vector opt.
254                  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)
255
256                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
257                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
258
259                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
260
261                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
262                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
263                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
264               END DO
265            END DO
266         ELSE                   ! other coordinate system (zco or sco) : e3t
267            DO jj = 2, jpj
268               DO ji = 1, fs_jpim1   ! vector opt.
269                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj)
270
271                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
272                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
273
274                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
275
276                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
277                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
278                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
279               END DO
280            END DO
281         ENDIF
282
283
284         ! Second derivative (divergence) and add to the general trend
285         ! -----------------------------------------------------------
286         DO jj = 2, jpjm1
287            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
288               ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
289                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
290               va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
291                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
292            END DO
293         END DO
294         !                                             ! ===============
295      END DO                                           !   End of slab
296      !                                                ! ===============
297
298      ! print sum trends (used for debugging)
299      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
300         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
301
302
303      !                                                ! ===============
304      DO jj = 2, jpjm1                                 !  Vertical slab
305         !                                             ! ===============
306
307 
308         ! I. vertical trends associated with the lateral mixing
309         ! =====================================================
310         !  (excluding the vertical flux proportional to dk[t]
311
312
313         ! I.1 horizontal momentum gradient
314         ! --------------------------------
315
316         DO jk = 1, jpk
317            DO ji = 2, jpi
318               ! i-gradient of u at jj
319               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
320               ! j-gradient of u and v at jj
321               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
322               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
323               ! j-gradient of u and v at jj+1
324               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
325               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
326            END DO
327         END DO
328         DO jk = 1, jpk
329            DO ji = 1, jpim1
330               ! i-gradient of v at jj
331               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
332            END DO
333         END DO
334
335
336         ! I.2 Vertical fluxes
337         ! -------------------
338
339         ! Surface and bottom vertical fluxes set to zero
340         DO ji = 1, jpi
341            zfuw(ji, 1 ) = 0.e0
342            zfvw(ji, 1 ) = 0.e0
343            zfuw(ji,jpk) = 0.e0
344            zfvw(ji,jpk) = 0.e0
345         END DO
346
347         ! interior (2=<jk=<jpk-1) on U field
348         DO jk = 2, jpkm1
349            DO ji = 2, jpim1
350               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk)
351               !
352               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
353               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
354               !
355               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
356                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
357               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   &
358                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. )
359
360               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
361               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
362               ! vertical flux on u field
363               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
364                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
365                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
366                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
367               ! update avmu (add isopycnal vertical coefficient to avmu)
368               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
369               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
370            END DO
371         END DO
372
373         ! interior (2=<jk=<jpk-1) on V field
374         DO jk = 2, jpkm1
375            DO ji = 2, jpim1
376               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk)
377
378               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
379               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
380
381               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
382                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
383               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
384                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
385
386               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
387               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
388               ! vertical flux on v field
389               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
390                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
391                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
392                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
393               ! update avmv (add isopycnal vertical coefficient to avmv)
394               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
395               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
396            END DO
397         END DO
398
399
400         ! I.3 Divergence of vertical fluxes added to the general tracer trend
401         ! -------------------------------------------------------------------
402         DO jk = 1, jpkm1
403            DO ji = 2, jpim1
404               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
405               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
406            END DO
407         END DO
408         !                                             ! ===============
409      END DO                                           !   End of slab
410      !                                                ! ===============
411      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
412      !
413      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso')
414      !
415   END SUBROUTINE dyn_ldf_iso
416
417   !!======================================================================
418END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.