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.
dynzdf.F90 in NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90 @ 10030

Last change on this file since 10030 was 10030, checked in by gm, 5 years ago

#1911 (ENHANCE-04): RK3 branch - step II.3 remove e3uw_$ e3vw_$, except e3.w_0 and use only after e3 in dyn/trazdf
  • Property svn:keywords set to Id
File size: 20.8 KB
Line 
1MODULE dynzdf
2   !!==============================================================================
3   !!                 ***  MODULE  dynzdf  ***
4   !! Ocean dynamics :  vertical component of the momentum mixing trend
5   !!==============================================================================
6   !! History :  1.0  !  2005-11  (G. Madec)  Original code
7   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
8   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
13   !!       zdf_trd   : diagnose the zdf velocity trends and the KE dissipation trend
14!!gm                        ==>>> zdf_trd currently not used
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE phycst         ! physical constants
18   USE dom_oce        ! ocean space and time domain variables
19   USE sbc_oce        ! surface boundary condition: ocean
20   USE zdf_oce        ! ocean vertical physics variables
21   USE zdfdrg         ! vertical physics: top/bottom drag coef.
22   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
23   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
24   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
25   USE trd_oce        ! trends: ocean variables
26   USE trddyn         ! trend manager: dynamics
27   !
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE iom             ! IOM library
31   USE prtctl         ! Print control
32   USE timing         ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   dyn_zdf   !  routine called by step.F90
38
39   !! * Substitutions
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
43   !! $Id$
44   !! Software governed by the CeCILL licence (./LICENSE)
45   !!----------------------------------------------------------------------
46CONTAINS
47   
48   SUBROUTINE dyn_zdf( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dyn_zdf  ***
51      !!
52      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
53      !!              together with the Leap-Frog time stepping using an
54      !!              implicit scheme.
55      !!
56      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
57      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
58      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
59      !!               - update the after velocity with the implicit vertical mixing.
60      !!      This requires to solver the following system:
61      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
62      !!      with the following surface/top/bottom boundary condition:
63      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
64      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
65      !!
66      !! ** Action :   (ua,va)   after velocity
67      !!---------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
69      !
70      INTEGER  ::   ji, jj, jk    ! dummy loop indices
71      INTEGER  ::   iku, ikv      ! local integers
72      REAL(wp) ::   zzwi, zDt_2   ! local scalars
73      REAL(wp) ::   zzws          !   -      -
74      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace
75      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      -
76      REAL(wp)                     ::   ze3uw_A , ze3uw_Ap1   ! local real
77      REAL(wp)                     ::   ze3vw_A , ze3vw_Ap1   ! local real
78      REAL(wp), DIMENSION(jpi,jpj) ::   zsshu_hA, zsshv_hA    ! 2D workspace
79      !!---------------------------------------------------------------------
80      !
81      IF( ln_timing )   CALL timing_start('dyn_zdf')
82      !
83      IF( kt == nit000 ) THEN       !* initialization
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
86         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
87      ENDIF
88      !
89      zDt_2 = rDt * 0.5_wp          !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2
90      !
91      !
92      !                             !* explicit top/bottom drag case
93      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
94      !
95      !
96      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
97         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
98         ztrdu(:,:,:) = ua(:,:,:)
99         ztrdv(:,:,:) = va(:,:,:)
100      ENDIF
101      !
102      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
103      !
104      !                    ! time stepping except vertical diffusion
105      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
106         DO jk = 1, jpkm1
107            ua(:,:,jk) = ( ub(:,:,jk) + rDt * ua(:,:,jk) ) * umask(:,:,jk)
108            va(:,:,jk) = ( vb(:,:,jk) + rDt * va(:,:,jk) ) * vmask(:,:,jk)
109         END DO
110      ELSE                                      ! applied on thickness weighted velocity
111         DO jk = 1, jpkm1
112            ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + rDt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
113            va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + rDt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
114         END DO
115      ENDIF
116      !                    ! add top/bottom friction
117      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
118      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
119      !                not lead to the effective stress seen over the whole barotropic loop.
120      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
121      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
122         DO jk = 1, jpkm1        ! remove barotropic velocities
123            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
124            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
125         END DO
126         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
127            DO ji = fs_2, fs_jpim1   ! vector opt.
128               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
129               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
130               ua(ji,jj,iku) = ua(ji,jj,iku) + zDt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / e3u_a(ji,jj,iku)
131               va(ji,jj,ikv) = va(ji,jj,ikv) + zDt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / e3v_a(ji,jj,ikv)
132            END DO
133         END DO
134         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
135            DO jj = 2, jpjm1       
136               DO ji = fs_2, fs_jpim1   ! vector opt.
137                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
138                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
139                  ua(ji,jj,iku) = ua(ji,jj,iku) + zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / e3u_a(ji,jj,iku)
140                  va(ji,jj,ikv) = va(ji,jj,ikv) + zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / e3v_a(ji,jj,ikv)
141               END DO
142            END DO
143         END IF
144      ENDIF
145      !
146      !              !==  Vertical diffusion on u  ==!
147      !
148      !                          !*  multiplicative factors on e3uw(Naa) and e3vw(Naa)
149      !
150      IF( ln_linssh ) THEN             !--  linear ssh case
151         DO jj = 1, jpjm1
152            DO ji = 1, jpim1
153               zsshu_hA(ji,jj) = 0._wp       ! no time variation in e3
154               zsshv_hA(ji,jj) = 0._wp
155            END DO
156         END DO
157      ELSE                             !--  Non linear ssh case
158         DO jj = 1, jpjm1
159            DO ji = 1, jpim1
160               zsshu_hA(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj  ,Naa) ) * r1_hu_0(ji,jj) * ssumask(ji,jj)
161               zsshv_hA(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji  ,jj+1,Naa) ) * r1_hv_0(ji,jj) * ssvmask(ji,jj)
162            END DO
163         END DO
164      ENDIF
165
166
167
168
169
170      SELECT CASE( nldf_dyn )    !* Matrix construction
171      !
172      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
173         DO jk = 1, jpkm1
174            DO jj = 2, jpjm1 
175               DO ji = fs_2, fs_jpim1   ! vector opt.
176!!gm Note that below, since ze3uw_A is used in a expression masked by wumask,
177!!      one can remove wumask from its expression   (same for ze3uw_Ap1
178                  ze3uw_A   = e3uw_0(ji,jj,jk  ) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk  ) )
179                  ze3uw_Ap1 = e3uw_0(ji,jj,jk+1) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk+1) )
180                  !
181                  zzwi = - rDt * ( 0.5 * (   avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) +   akzu(ji,jj,jk  ) )   &
182                     &                 / ( e3u_a(ji  ,jj,jk  ) * ze3uw_A         ) * wumask(ji,jj,jk  )
183                  zzws = - rDt * ( 0.5 * (   avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) +   akzu(ji,jj,jk+1) )   &
184                     &                 / ( e3u_a(ji  ,jj,jk  ) * ze3uw_Ap1       ) * wumask(ji,jj,jk+1)
185                  zwi(ji,jj,jk) = zzwi
186                  zws(ji,jj,jk) = zzws
187                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
188               END DO
189            END DO
190         END DO
191      CASE DEFAULT                  ! iso-level lateral mixing
192         DO jk = 1, jpkm1
193            DO jj = 2, jpjm1 
194               DO ji = fs_2, fs_jpim1   ! vector opt.
195!!gm Note that below, since ze3uw_A is used in a expression masked by wumask,
196!!      one can remove wumask from its expression   (same for ze3uw_Ap1
197                  ze3uw_A   = e3uw_0(ji,jj,jk  ) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk  ) )
198                  ze3uw_Ap1 = e3uw_0(ji,jj,jk+1) * ( 1._wp + zsshu_hA(ji,jj) * wumask(ji,jj,jk+1) )
199                  !
200                  zzwi = - zDt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( e3u_a(ji,jj,jk) * ze3uw_A   ) * wumask(ji,jj,jk  )
201                  zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( e3u_a(ji,jj,jk) * ze3uw_Ap1 ) * wumask(ji,jj,jk+1)
202                  zwi(ji,jj,jk) = zzwi
203                  zws(ji,jj,jk) = zzws
204                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
205               END DO
206            END DO
207         END DO
208      END SELECT
209      !
210      DO jj = 2, jpjm1           !* Surface boundary conditions
211         DO ji = fs_2, fs_jpim1     ! vector opt.
212            zwi(ji,jj,1) = 0._wp
213            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
214         END DO
215      END DO
216      !
217      !              !==  Apply semi-implicit bottom friction  ==!
218      !
219      !     Only needed for semi-implicit bottom friction setup. The explicit
220      !     bottom friction has been included in "u(v)a" which act as the R.H.S
221      !     column vector of the tri-diagonal matrix equation
222      !
223      IF ( ln_drgimp ) THEN      ! implicit bottom friction
224         DO jj = 2, jpjm1
225            DO ji = 2, jpim1
226               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
227               zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_a(ji,jj,iku)
228            END DO
229         END DO
230         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
231            DO jj = 2, jpjm1
232               DO ji = 2, jpim1
233                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
234                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
235                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_a(ji,jj,iku)
236               END DO
237            END DO
238         END IF
239      ENDIF
240      !
241      ! Matrix inversion starting from the first level
242      !-----------------------------------------------------------------------
243      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
244      !
245      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
246      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
247      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
248      !        (        ...               )( ...  ) ( ...  )
249      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
250      !
251      !   m is decomposed in the product of an upper and a lower triangular matrix
252      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
253      !   The solution (the after velocity) is in ua
254      !-----------------------------------------------------------------------
255      !
256      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
257         DO jj = 2, jpjm1   
258            DO ji = fs_2, fs_jpim1   ! vector opt.
259               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
260            END DO
261         END DO
262      END DO
263      !
264      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
265         DO ji = fs_2, fs_jpim1   ! vector opt.
266            ua(ji,jj,1) = ua(ji,jj,1) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u_a(ji,jj,1) * rho0 ) * umask(ji,jj,1)
267         END DO
268      END DO
269      DO jk = 2, jpkm1
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1
272               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
273            END DO
274         END DO
275      END DO
276      !
277      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
278         DO ji = fs_2, fs_jpim1   ! vector opt.
279            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
280         END DO
281      END DO
282      DO jk = jpk-2, 1, -1
283         DO jj = 2, jpjm1
284            DO ji = fs_2, fs_jpim1
285               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
286            END DO
287         END DO
288      END DO
289      !
290      !              !==  Vertical diffusion on v  ==!
291      !
292      !                     
293      SELECT CASE( nldf_dyn )    !* Matrix construction
294      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
295         DO jk = 1, jpkm1
296            DO jj = 2, jpjm1   
297               DO ji = fs_2, fs_jpim1   ! vector opt.
298                  ze3vw_A   = e3vw_0(ji,jj,jk  ) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk  ) )
299                  ze3vw_Ap1 = e3vw_0(ji,jj,jk+1) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk+1) )
300                  !
301                  zzwi = - rDt * ( 0.5 * (   avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) +   akzv(ji,jj,jk  ) )   &
302                     &                 / ( e3v_a(ji,jj  ,jk  ) * ze3vw_A         ) * wvmask(ji,jj,jk  )
303                  zzws = - rDt * ( 0.5 * (   avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) +   akzv(ji,jj,jk+1) )   &
304                     &                 / ( e3v_a(ji,jj  ,jk  ) * ze3vw_Ap1       ) * wvmask(ji,jj,jk+1)
305                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
306                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
307                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
308               END DO
309            END DO
310         END DO
311      CASE DEFAULT                  ! iso-level lateral mixing
312         DO jk = 1, jpkm1
313            DO jj = 2, jpjm1   
314               DO ji = fs_2, fs_jpim1   ! vector opt.
315                  ze3vw_A   = e3vw_0(ji,jj,jk  ) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk  ) )
316                  ze3vw_Ap1 = e3vw_0(ji,jj,jk+1) * ( 1._wp + zsshv_hA(ji,jj) * wvmask(ji,jj,jk+1) )
317                  !
318                  zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( e3v_a(ji,jj,jk) * ze3vw_A   ) * wvmask(ji,jj,jk  )
319                  zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( e3v_a(ji,jj,jk) * ze3vw_Ap1 ) * wvmask(ji,jj,jk+1)
320                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
321                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
322                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
323               END DO
324            END DO
325         END DO
326      END SELECT
327      !
328      DO jj = 2, jpjm1           !* Surface boundary conditions
329         DO ji = fs_2, fs_jpim1     ! vector opt.
330            zwi(ji,jj,1) = 0._wp
331            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
332         END DO
333      END DO
334      !              !==  Apply semi-implicit top/bottom friction  ==!
335      !
336      !     Only needed for semi-implicit bottom friction setup. The explicit
337      !     bottom friction has been included in "u(v)a" which act as the R.H.S
338      !     column vector of the tri-diagonal matrix equation
339      !
340      IF( ln_drgimp ) THEN
341         DO jj = 2, jpjm1
342            DO ji = 2, jpim1
343               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
344               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_a(ji,jj,ikv)           
345            END DO
346         END DO
347         IF ( ln_isfcav ) THEN
348            DO jj = 2, jpjm1
349               DO ji = 2, jpim1
350                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
351                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3v_a(ji,jj,ikv)
352               END DO
353            END DO
354         ENDIF
355      ENDIF
356
357      ! Matrix inversion
358      !-----------------------------------------------------------------------
359      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
360      !
361      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
362      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
363      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
364      !        (        ...               )( ...  ) ( ...  )
365      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
366      !
367      !   m is decomposed in the product of an upper and lower triangular matrix
368      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
369      !   The solution (after velocity) is in 2d array va
370      !-----------------------------------------------------------------------
371      !
372      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
373         DO jj = 2, jpjm1   
374            DO ji = fs_2, fs_jpim1   ! vector opt.
375               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
376            END DO
377         END DO
378      END DO
379      !
380      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
381         DO ji = fs_2, fs_jpim1   ! vector opt.         
382            va(ji,jj,1) = va(ji,jj,1) + zDt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v_a(ji,jj,1) * rho0 ) * vmask(ji,jj,1)
383         END DO
384      END DO
385      DO jk = 2, jpkm1
386         DO jj = 2, jpjm1
387            DO ji = fs_2, fs_jpim1   ! vector opt.
388               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
389            END DO
390         END DO
391      END DO
392      !
393      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
394         DO ji = fs_2, fs_jpim1   ! vector opt.
395            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
396         END DO
397      END DO
398      DO jk = jpk-2, 1, -1
399         DO jj = 2, jpjm1
400            DO ji = fs_2, fs_jpim1
401               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
402            END DO
403         END DO
404      END DO
405      !
406      IF( l_trddyn ) THEN                        ! save the vertical diffusive trends for further diagnostics
407         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
408            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ztrdu(:,:,:)
409            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ztrdv(:,:,:)
410         ELSE                                      ! applied on thickness weighted velocity
411            ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ztrdu(:,:,:)
412            ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ztrdv(:,:,:)
413         ENDIF
414         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
415         DEALLOCATE( ztrdu, ztrdv ) 
416      ENDIF
417      !                                          ! print mean trends (used for debugging)
418      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
419         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
420         !
421      IF( ln_timing )   CALL timing_stop('dyn_zdf')
422      !
423   END SUBROUTINE dyn_zdf
424
425   !!==============================================================================
426END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.