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/RK3_SRC/DYN – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DYN/dynldf_iso.F90 @ 8367

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

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