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

source: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90 @ 13226

Last change on this file since 13226 was 13226, checked in by orioltp, 4 years ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

  • Property svn:keywords set to Id
File size: 20.1 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 "do_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_3D_00_00( 1, jpk )
130            uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
131            vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
132            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
133            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
134         END_3D
135         ! Lateral boundary conditions on the slopes
136         CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
137         !
138       ENDIF
139         
140      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
141     
142      !                                                ! ===============
143      DO jk = 1, jpkm1                                 ! Horizontal slab
144         !                                             ! ===============
145
146         ! Vertical u- and v-shears at level jk and jk+1
147         ! ---------------------------------------------
148         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
149         !                             zdkv(jk=1)=zdkv(jk=2)
150
151         zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
152         zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
153
154         IF( jk == 1 ) THEN
155            zdku(:,:) = zdk1u(:,:)
156            zdkv(:,:) = zdk1v(:,:)
157         ELSE
158            zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
159            zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
160         ENDIF
161
162         !                               -----f-----
163         ! Horizontal fluxes on U             | 
164         ! --------------------===        t   u   t
165         !                                    | 
166         ! i-flux at t-point             -----f-----
167
168         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
169            DO_2D_00_01
170               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)
171
172               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
173                  &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
174
175               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
176
177               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
178                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
179                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
180            END_2D
181         ELSE                   ! other coordinate system (zco or sco) : e3t
182            DO_2D_00_01
183               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
184
185               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
186                  &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
187
188               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
189
190               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
191                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
192                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
193            END_2D
194         ENDIF
195
196         ! j-flux at f-point
197         DO_2D_10_10
198            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
199
200            zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
201               &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
202
203            zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
204
205            zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
206               &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
207               &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
208         END_2D
209
210         !                                |   t   |
211         ! Horizontal fluxes on V         |       |
212         ! --------------------===        f---v---f
213         !                                |       |
214         ! i-flux at f-point              |   t   |
215
216         DO_2D_00_10
217            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
218
219            zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
220               &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
221
222            zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
223
224            zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
225               &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
226               &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
227         END_2D
228
229         ! j-flux at t-point
230         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
231            DO_2D_01_10
232               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)
233
234               zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
235                  &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
236
237               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
238
239               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
240                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
241                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
242            END_2D
243         ELSE                   ! other coordinate system (zco or sco) : e3t
244            DO_2D_01_10
245               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
246
247               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
248                  &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
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_2D
256         ENDIF
257
258
259         ! Second derivative (divergence) and add to the general trend
260         ! -----------------------------------------------------------
261         DO_2D_00_00
262            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
263               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
264            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
265               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
266         END_2D
267         !                                             ! ===============
268      END DO                                           !   End of slab
269      !                                                ! ===============
270
271      ! print sum trends (used for debugging)
272      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
273         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
274
275
276      !                                                ! ===============
277      DO jj = 2, jpjm1                                 !  Vertical slab
278         !                                             ! ===============
279
280 
281         ! I. vertical trends associated with the lateral mixing
282         ! =====================================================
283         !  (excluding the vertical flux proportional to dk[t]
284
285
286         ! I.1 horizontal momentum gradient
287         ! --------------------------------
288
289         DO jk = 1, jpk
290            DO ji = 2, jpi
291               ! i-gradient of u at jj
292               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) )
293               ! j-gradient of u and v at jj
294               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) )
295               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) )
296               ! j-gradient of u and v at jj+1
297               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) )
298               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
299            END DO
300         END DO
301         DO jk = 1, jpk
302            DO ji = 1, jpim1
303               ! i-gradient of v at jj
304               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
305            END DO
306         END DO
307
308
309         ! I.2 Vertical fluxes
310         ! -------------------
311
312         ! Surface and bottom vertical fluxes set to zero
313         DO ji = 1, jpi
314            zfuw(ji, 1 ) = 0.e0
315            zfvw(ji, 1 ) = 0.e0
316            zfuw(ji,jpk) = 0.e0
317            zfvw(ji,jpk) = 0.e0
318         END DO
319
320         ! interior (2=<jk=<jpk-1) on U field
321         DO jk = 2, jpkm1
322            DO ji = 2, jpim1
323               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
324               !
325               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
326               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
327               !
328               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
329                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
330               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
331                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
332
333               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
334               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
335               ! vertical flux on u field
336               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
337                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
338                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
339                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
340               ! vertical mixing coefficient (akzu)
341               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
342               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
343            END DO
344         END DO
345
346         ! interior (2=<jk=<jpk-1) on V field
347         DO jk = 2, jpkm1
348            DO ji = 2, jpim1
349               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
350               !
351               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
352               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
353               !
354               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
355                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
356               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
357                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
358
359               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
360               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
361               ! vertical flux on v field
362               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
363                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
364                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
365                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
366               ! vertical mixing coefficient (akzv)
367               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
368               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
369            END DO
370         END DO
371
372
373         ! I.3 Divergence of vertical fluxes added to the general tracer trend
374         ! -------------------------------------------------------------------
375         DO jk = 1, jpkm1
376            DO ji = 2, jpim1
377               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)
378               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)
379            END DO
380         END DO
381         !                                             ! ===============
382      END DO                                           !   End of slab
383      !                                                ! ===============
384   END SUBROUTINE dyn_ldf_iso
385
386   !!======================================================================
387END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.