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/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN/dynldf_iso.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago

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