source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OCE/DYN/dynldf_iso.F90 @ 11872

Last change on this file since 11872 was 11872, checked in by acc, 11 months ago

Branch 2019/fix_sn_cfctl_ticket2328. See #2328. Replacement of ln_ctl and activation of full functionality with
sn_cfctl structure. These changes rename structure components l_mppout and l_mpptop as l_prtctl and l_prttrc
and introduce l_glochk to activate former ln_ctl code in stpctl.F90 to perform global location of min and max
checks. Also added is l_allon which can be used to activate all output (much like the former ln_ctl). If l_allon
is .false. then l_config decides whether or not the suboptions are used.

   sn_cfctl%l_glochk = .FALSE.    ! Range sanity checks are local (F) or global (T). Set T for debugging only
   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
   sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the remaining options

Note, these changes pass SETTE tests but all references to ln_ctl need to be removed from the sette scripts.

  • Property svn:keywords set to Id
File size: 20.3 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 )
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      !!       -(ua,va) 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      !
107      INTEGER  ::   ji, jj, jk   ! dummy loop indices
108      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
109      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
110      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      -
111      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
112      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
113      !!----------------------------------------------------------------------
114      !
115      IF( kt == nit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
119         !                                      ! allocate dyn_ldf_bilap arrays
120         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
121      ENDIF
122
123!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
124      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
125      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
126         !
127         DO jk = 1, jpk         ! set the slopes of iso-level
128            DO jj = 2, jpjm1
129               DO ji = 2, jpim1
130                  uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
131                  vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
132                  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
133                  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
134               END DO
135            END DO
136         END DO
137         ! Lateral boundary conditions on the slopes
138         CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1., vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1. )
139         !
140       ENDIF
141         
142      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
143     
144      !                                                ! ===============
145      DO jk = 1, jpkm1                                 ! Horizontal slab
146         !                                             ! ===============
147
148         ! Vertical u- and v-shears at level jk and jk+1
149         ! ---------------------------------------------
150         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
151         !                             zdkv(jk=1)=zdkv(jk=2)
152
153         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
154         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
155
156         IF( jk == 1 ) THEN
157            zdku(:,:) = zdk1u(:,:)
158            zdkv(:,:) = zdk1v(:,:)
159         ELSE
160            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
161            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
162         ENDIF
163
164         !                               -----f-----
165         ! Horizontal fluxes on U             | 
166         ! --------------------===        t   u   t
167         !                                    | 
168         ! i-flux at t-point             -----f-----
169
170         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
171            DO jj = 2, jpjm1
172               DO ji = fs_2, jpi   ! vector opt.
173                  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)
174
175                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
176                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
177
178                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
179   
180                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
181                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
182                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
183               END DO
184            END DO
185         ELSE                   ! other coordinate system (zco or sco) : e3t
186            DO jj = 2, jpjm1
187               DO ji = fs_2, jpi   ! vector opt.
188                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj)
189
190                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
191                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
192
193                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
194
195                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
196                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
197                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
198               END DO
199            END DO
200         ENDIF
201
202         ! j-flux at f-point
203         DO jj = 1, jpjm1
204            DO ji = 1, fs_jpim1   ! vector opt.
205               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj)
206
207               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
208                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
209
210               zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
211
212               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
213                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
214                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
215            END DO
216         END DO
217
218         !                                |   t   |
219         ! Horizontal fluxes on V         |       |
220         ! --------------------===        f---v---f
221         !                                |       |
222         ! i-flux at f-point              |   t   |
223
224         DO jj = 2, jpjm1
225            DO ji = 1, fs_jpim1   ! vector opt.
226               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj)
227
228               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
229                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
230
231               zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
232
233               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
234                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
235                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
236            END DO
237         END DO
238
239         ! j-flux at t-point
240         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
241            DO jj = 2, jpj
242               DO ji = 1, fs_jpim1   ! vector opt.
243                  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)
244
245                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
246                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
247
248                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
249
250                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
251                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
252                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
253               END DO
254            END DO
255         ELSE                   ! other coordinate system (zco or sco) : e3t
256            DO jj = 2, jpj
257               DO ji = 1, fs_jpim1   ! vector opt.
258                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj)
259
260                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
261                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
262
263                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
264
265                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
266                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
267                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
268               END DO
269            END DO
270         ENDIF
271
272
273         ! Second derivative (divergence) and add to the general trend
274         ! -----------------------------------------------------------
275         DO jj = 2, jpjm1
276            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
277               ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
278                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
279               va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
280                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
281            END DO
282         END DO
283         !                                             ! ===============
284      END DO                                           !   End of slab
285      !                                                ! ===============
286
287      ! print sum trends (used for debugging)
288      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
289         &                                  tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
290
291
292      !                                                ! ===============
293      DO jj = 2, jpjm1                                 !  Vertical slab
294         !                                             ! ===============
295
296 
297         ! I. vertical trends associated with the lateral mixing
298         ! =====================================================
299         !  (excluding the vertical flux proportional to dk[t]
300
301
302         ! I.1 horizontal momentum gradient
303         ! --------------------------------
304
305         DO jk = 1, jpk
306            DO ji = 2, jpi
307               ! i-gradient of u at jj
308               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
309               ! j-gradient of u and v at jj
310               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
311               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
312               ! j-gradient of u and v at jj+1
313               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
314               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
315            END DO
316         END DO
317         DO jk = 1, jpk
318            DO ji = 1, jpim1
319               ! i-gradient of v at jj
320               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
321            END DO
322         END DO
323
324
325         ! I.2 Vertical fluxes
326         ! -------------------
327
328         ! Surface and bottom vertical fluxes set to zero
329         DO ji = 1, jpi
330            zfuw(ji, 1 ) = 0.e0
331            zfvw(ji, 1 ) = 0.e0
332            zfuw(ji,jpk) = 0.e0
333            zfvw(ji,jpk) = 0.e0
334         END DO
335
336         ! interior (2=<jk=<jpk-1) on U field
337         DO jk = 2, jpkm1
338            DO ji = 2, jpim1
339               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
340               !
341               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
342               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
343               !
344               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
345                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
346               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
347                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
348
349               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
350               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
351               ! vertical flux on u field
352               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
353                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
354                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
355                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
356               ! vertical mixing coefficient (akzu)
357               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
358               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
359            END DO
360         END DO
361
362         ! interior (2=<jk=<jpk-1) on V field
363         DO jk = 2, jpkm1
364            DO ji = 2, jpim1
365               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
366               !
367               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
368               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
369               !
370               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
371                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
372               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
373                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
374
375               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
376               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
377               ! vertical flux on v field
378               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
379                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
380                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
381                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
382               ! vertical mixing coefficient (akzv)
383               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
384               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
385            END DO
386         END DO
387
388
389         ! I.3 Divergence of vertical fluxes added to the general tracer trend
390         ! -------------------------------------------------------------------
391         DO jk = 1, jpkm1
392            DO ji = 2, jpim1
393               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
394               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
395            END DO
396         END DO
397         !                                             ! ===============
398      END DO                                           !   End of slab
399      !                                                ! ===============
400   END SUBROUTINE dyn_ldf_iso
401
402   !!======================================================================
403END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.