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_lf.F90 in NEMO/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynldf_iso_lf.F90

Last change on this file was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling
File size: 21.4 KB
Line 
1MODULE dynldf_iso_lf
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_lf           ! called by step.F90
35   PUBLIC   dyn_ldf_iso_alloc_lf     ! called by nemogcm.F90
36
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id: dynldf_iso.F90 14757 2021-04-27 15:33:44Z francesca $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   INTEGER FUNCTION dyn_ldf_iso_alloc_lf()
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
52      !!----------------------------------------------------------------------
53      dyn_ldf_iso_alloc_lf = 0
54      IF( .NOT. ALLOCATED( akzu ) ) THEN
55         ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,jpk), STAT=dyn_ldf_iso_alloc_lf )
56            !
57         IF( dyn_ldf_iso_alloc_lf /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
58      ENDIF
59   END FUNCTION dyn_ldf_iso_alloc_lf
60
61
62   SUBROUTINE dyn_ldf_iso_lf( 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) ::   zdiu, zdiu_km1, zdiu_ip1, zdiu_ip1_km1       !   -      -
114      REAL(wp) ::   zdju, zdju_km1, zdj1u, zdj1u_km1
115      REAL(wp) ::   zdjv, zdjv_km1, zdj1v, zdj1v_km1
116      REAL(wp) ::   zdiv_im1_km1, zdiv, zdiv_im1, zdiv_km1       !   -      -
117      REAL(wp), DIMENSION(A2D(nn_hls)) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
118      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
119      REAL(wp), DIMENSION(A1Di(nn_hls),jpk) ::   zfuw, zfvw
120      !!----------------------------------------------------------------------
121      !
122      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
123         IF( kt == nit000 ) THEN
124            IF(lwp) WRITE(numout,*)
125            IF(lwp) WRITE(numout,*) 'dyn_ldf_iso_lf : iso-neutral laplacian diffusive operator or '
126            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~     s-coordinate horizontal diffusive operator'
127            !                                      ! allocate dyn_ldf_bilap arrays
128            IF( dyn_ldf_iso_alloc_lf() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
129         ENDIF
130      ENDIF
131
132!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
133      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
134      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
135         !
136         DO_3D_OVR( 1, 1, 1, 1, 1, jpk )      ! set the slopes of iso-level
137            uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
138            vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
139            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
140            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
141         END_3D
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         DO_2D( 1, 1, 1, 1 )
157            zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1)
158            zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1)
159         END_2D
160
161         IF( jk == 1 ) THEN
162            zdku(:,:) = zdk1u(:,:)
163            zdkv(:,:) = zdk1v(:,:)
164         ELSE
165            DO_2D( 1, 1, 1, 1 )
166               zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
167               zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
168            END_2D
169         ENDIF
170
171         !                               -----f-----
172         ! Horizontal fluxes on U             | 
173         ! --------------------===        t   u   t
174         !                                    | 
175         ! i-flux at t-point             -----f-----
176
177         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
178            DO_2D( 0, 1, 0, 0 )
179               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
180                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
181                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
182
183               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
184                  &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
185
186               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
187
188               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
189                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
190                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
191            END_2D
192         ELSE                   ! other coordinate system (zco or sco) : e3t
193            DO_2D( 0, 1, 0, 0 )
194               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
195                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
196
197               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
198                  &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
199
200               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
201
202               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
203                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
204                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
205            END_2D
206         ENDIF
207
208         ! j-flux at f-point
209         DO_2D( 1, 0, 1, 0 )
210            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
211               &     * e1f(ji,jj) * e3f(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 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
217
218            zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
219               &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
220               &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
221
222            !                                |   t   |
223            ! Horizontal fluxes on V         |       |
224            ! --------------------===        f---v---f
225            !                                |       |
226            ! i-flux at f-point              |   t   |
227
228            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
229               &     * e2f(ji,jj) * e3f(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 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
237               &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
238               &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
239         END_2D
240
241         ! j-flux at t-point
242         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
243            DO_2D( 1, 0, 0, 1 )
244               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
245                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
246                  &            e3v(ji,jj-1,jk,Kmm) ) * 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 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
254                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
255                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
256            END_2D
257         ELSE                   ! other coordinate system (zco or sco) : e3t
258            DO_2D( 1, 0, 0, 1 )
259               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
260                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
261
262               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
263                  &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
264
265               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
266
267               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
268                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
269                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
270            END_2D
271         ENDIF
272
273
274         ! Second derivative (divergence) and add to the general trend
275         ! -----------------------------------------------------------
276         DO_2D( 0, 0, 0, 0 )      !!gm Question vectop possible??? !!bug
277            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
278               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
279               &                           / e3u(ji,jj,jk,Kmm)
280            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
281               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   &
282               &                           / e3v(ji,jj,jk,Kmm)
283         END_2D
284         !                                             ! ===============
285      END DO                                           !   End of slab
286      !                                                ! ===============
287
288      ! print sum trends (used for debugging)
289      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
290         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
291
292
293      !                                                ! ===============
294      DO jj = ntsj, ntej                               !  Vertical slab
295         !                                             ! ===============
296
297 
298         ! I. vertical trends associated with the lateral mixing
299         ! =====================================================
300         !  (excluding the vertical flux proportional to dk[t]
301
302         ! I.2 Vertical fluxes
303         ! -------------------
304
305         ! Surface and bottom vertical fluxes set to zero
306         DO ji = ntsi - nn_hls, ntei + nn_hls
307            zfuw(ji, 1 ) = 0.e0
308            zfvw(ji, 1 ) = 0.e0
309            zfuw(ji,jpk) = 0.e0
310            zfvw(ji,jpk) = 0.e0
311         END DO
312
313         ! interior (2=<jk=<jpk-1) on U and V fields
314         DO jk = 2, jpkm1
315            DO ji = ntsi, ntei
316               ! I.1 horizontal momentum gradient
317               ! --------------------------------
318               ! i-gradient of u at jj
319               zdiu = tmask(ji,jj,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) )
320               zdiu_km1 = tmask(ji,jj,jk-1) * ( puu(ji,jj,jk-1,Kbb) - puu(ji-1,jj,jk-1,Kbb) )
321               zdiu_ip1 = tmask(ji+1,jj,jk) * ( puu(ji+1,jj,jk,Kbb) - puu(ji,jj,jk,Kbb) )
322               zdiu_ip1_km1 = tmask(ji+1,jj,jk-1) * ( puu(ji+1,jj,jk-1,Kbb) - puu(ji,jj,jk-1,Kbb) )
323               ! j-gradient of u and v at jj
324               zdju = fmask(ji,jj,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )
325               zdju_km1 = fmask(ji,jj,jk-1) * ( puu(ji,jj+1,jk-1,Kbb) - puu(ji,jj,jk-1,Kbb) )
326               ! j-gradient of u and v at jj+1
327               zdj1u = fmask(ji,jj-1,jk) * ( puu(ji,jj,jk,Kbb) - puu(ji,jj-1,jk,Kbb) )
328               zdj1u_km1 = fmask(ji,jj-1,jk-1) * ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj-1,jk-1,Kbb) )
329               !
330               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
331               !
332               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
333               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
334               !
335               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
336                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
337               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
338                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
339
340               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
341               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
342               ! vertical flux on u field
343               zfuw(ji,jk) = zcof3 * (  zdiu_km1 + zdiu_ip1_km1 + zdiu + zdiu_ip1 )   &
344                  &        + zcof4 * (  zdj1u_km1 + zdju_km1 + zdj1u + zdju )
345               ! vertical mixing coefficient (akzu)
346               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
347               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
348
349               ! I.1 horizontal momentum gradient
350               ! --------------------------------
351               ! j-gradient of u and v at jj
352               zdjv     = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) )
353               zdjv_km1 = tmask(ji,jj,jk-1) * ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj-1,jk-1,Kbb) )
354               ! i-gradient of v at jj
355               zdiv = fmask(ji,jj,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )
356               zdiv_im1 = fmask(ji-1,jj,jk) * ( pvv(ji,jj,jk,Kbb) - pvv(ji-1,jj,jk,Kbb) )
357               zdiv_km1 = fmask(ji,jj,jk-1) * ( pvv(ji+1,jj,jk-1,Kbb) - pvv(ji,jj,jk-1,Kbb) )
358               zdiv_im1_km1 = fmask(ji-1,jj,jk-1) * ( pvv(ji,jj,jk-1,Kbb) - pvv(ji-1,jj,jk-1,Kbb) )
359               ! j-gradient of u and v at jj+1
360               zdj1v = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji,jj,jk,Kbb) )
361               zdj1v_km1 = tmask(ji,jj+1,jk-1) * ( pvv(ji,jj+1,jk-1,Kbb) - pvv(ji,jj,jk-1,Kbb) )
362               !
363               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
364               !
365               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
366               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
367               !
368               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
369                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
370               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
371                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
372
373               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
374               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
375               ! vertical flux on v field
376               zfvw(ji,jk) = zcof3 * (  zdiv_km1 + zdiv_im1_km1 + zdiv + zdiv_im1 )   &
377                  &        + zcof4 * (  zdjv_km1 + zdj1v_km1 + zdjv + zdj1v )
378               ! vertical mixing coefficient (akzv)
379               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
380               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
381            END DO
382         END DO
383
384
385         ! I.3 Divergence of vertical fluxes added to the general tracer trend
386         ! -------------------------------------------------------------------
387         DO jk = 1, jpkm1
388            DO ji = ntsi, ntei
389               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   &
390                  &               / e3u(ji,jj,jk,Kmm)
391               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   &
392                  &               / e3v(ji,jj,jk,Kmm)
393            END DO
394         END DO
395         !                                             ! ===============
396      END DO                                           !   End of slab
397      !                                                ! ===============
398   END SUBROUTINE dyn_ldf_iso_lf
399
400   !!======================================================================
401END MODULE dynldf_iso_lf
Note: See TracBrowser for help on using the repository browser.