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

source: NEMO/trunk/src/OCE/DYN/dynzdf.F90 @ 10364

Last change on this file since 10364 was 10364, checked in by acc, 5 years ago

Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.
  • Property svn:keywords set to Id
File size: 26.1 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   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE phycst         ! physical constants
16   USE dom_oce        ! ocean space and time domain variables
17   USE sbc_oce        ! surface boundary condition: ocean
18   USE zdf_oce        ! ocean vertical physics variables
19   USE zdfdrg         ! vertical physics: top/bottom drag coef.
20   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
21   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
22   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
23   USE trd_oce        ! trends: ocean variables
24   USE trddyn         ! trend manager: dynamics
25   !
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! MPP library
28   USE prtctl         ! Print control
29   USE timing         ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_zdf   !  routine called by step.F90
35
36   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46   
47   SUBROUTINE dyn_zdf( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_zdf  ***
50      !!
51      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
52      !!              together with the Leap-Frog time stepping using an
53      !!              implicit scheme.
54      !!
55      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
56      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
57      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
58      !!               - update the after velocity with the implicit vertical mixing.
59      !!      This requires to solver the following system:
60      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
61      !!      with the following surface/top/bottom boundary condition:
62      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
63      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
64      !!
65      !! ** Action :   (ua,va)   after velocity
66      !!---------------------------------------------------------------------
67      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
68      !
69      INTEGER  ::   ji, jj, jk         ! dummy loop indices
70      INTEGER  ::   iku, ikv           ! local integers
71      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
72      REAL(wp) ::   zzws, ze3va        !   -      -
73      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
74      REAL(wp) ::   zWu , zWv          !   -      -
75      REAL(wp) ::   zWui, zWvi         !   -      -
76      REAL(wp) ::   zWus, zWvs         !   -      -
77      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
78      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
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         !
88         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
89         ELSE                   ;    r_vvl = 1._wp
90         ENDIF
91      ENDIF
92      !                             !* set time step
93      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
94      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
95      ENDIF
96      !
97      !                             !* explicit top/bottom drag case
98      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
99      !
100      !
101      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
102         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
103         ztrdu(:,:,:) = ua(:,:,:)
104         ztrdv(:,:,:) = va(:,:,:)
105      ENDIF
106      !
107      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
108      !
109      !                    ! time stepping except vertical diffusion
110      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
111         DO jk = 1, jpkm1
112            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
113            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
114         END DO
115      ELSE                                      ! applied on thickness weighted velocity
116         DO jk = 1, jpkm1
117            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
118               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
119            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
120               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
121         END DO
122      ENDIF
123      !                    ! add top/bottom friction
124      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
125      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
126      !                not lead to the effective stress seen over the whole barotropic loop.
127      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
128      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
129         DO jk = 1, jpkm1        ! remove barotropic velocities
130            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
131            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
132         END DO
133         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
134            DO ji = fs_2, fs_jpim1   ! vector opt.
135               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
136               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
137               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
138               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
139               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
140               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
141            END DO
142         END DO
143         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
144            DO jj = 2, jpjm1       
145               DO ji = fs_2, fs_jpim1   ! vector opt.
146                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
147                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
148                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
149                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
150                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
151                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
152               END DO
153            END DO
154         END IF
155      ENDIF
156      !
157      !              !==  Vertical diffusion on u  ==!
158      !
159      !                    !* Matrix construction
160      zdt = r2dt * 0.5
161      IF( ln_zad_Aimp ) THEN   !!
162         SELECT CASE( nldf_dyn )
163         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
164            DO jk = 1, jpkm1
165               DO jj = 2, jpjm1 
166                  DO ji = fs_2, fs_jpim1   ! vector opt.
167                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
168                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
169                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
170                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
171                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
172                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
173                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
174                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
175                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
176                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
177                  END DO
178               END DO
179            END DO
180         CASE DEFAULT               ! iso-level lateral mixing
181            DO jk = 1, jpkm1
182               DO jj = 2, jpjm1 
183                  DO ji = fs_2, fs_jpim1   ! vector opt.
184                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
185                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
186                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
187                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
188                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
189                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
190                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
191                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
192                  END DO
193               END DO
194            END DO
195         END SELECT
196         DO jj = 2, jpjm1     !* Surface boundary conditions
197            DO ji = fs_2, fs_jpim1   ! vector opt.
198               zwi(ji,jj,1) = 0._wp
199               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
200               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
201               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) )
202               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
203               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
204            END DO
205         END DO
206      ELSE
207         SELECT CASE( nldf_dyn )
208         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
209            DO jk = 1, jpkm1
210               DO jj = 2, jpjm1 
211                  DO ji = fs_2, fs_jpim1   ! vector opt.
212                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
213                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
214                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
215                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
216                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
217                     zwi(ji,jj,jk) = zzwi
218                     zws(ji,jj,jk) = zzws
219                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
220                  END DO
221               END DO
222            END DO
223         CASE DEFAULT               ! iso-level lateral mixing
224            DO jk = 1, jpkm1
225               DO jj = 2, jpjm1 
226                  DO ji = fs_2, fs_jpim1   ! vector opt.
227                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
228                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
229                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
230                     zwi(ji,jj,jk) = zzwi
231                     zws(ji,jj,jk) = zzws
232                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
233                  END DO
234               END DO
235            END DO
236         END SELECT
237         DO jj = 2, jpjm1     !* Surface boundary conditions
238            DO ji = fs_2, fs_jpim1   ! vector opt.
239               zwi(ji,jj,1) = 0._wp
240               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
241            END DO
242         END DO
243      ENDIF
244      !
245      !
246      !              !==  Apply semi-implicit bottom friction  ==!
247      !
248      !     Only needed for semi-implicit bottom friction setup. The explicit
249      !     bottom friction has been included in "u(v)a" which act as the R.H.S
250      !     column vector of the tri-diagonal matrix equation
251      !
252      IF ( ln_drgimp ) THEN      ! implicit bottom friction
253         DO jj = 2, jpjm1
254            DO ji = 2, jpim1
255               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
256               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
257               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
258            END DO
259         END DO
260         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
261            DO jj = 2, jpjm1
262               DO ji = 2, jpim1
263                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
264                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
265                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
266                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
267               END DO
268            END DO
269         END IF
270      ENDIF
271      !
272      ! Matrix inversion starting from the first level
273      !-----------------------------------------------------------------------
274      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
275      !
276      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
277      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
278      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
279      !        (        ...               )( ...  ) ( ...  )
280      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
281      !
282      !   m is decomposed in the product of an upper and a lower triangular matrix
283      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
284      !   The solution (the after velocity) is in ua
285      !-----------------------------------------------------------------------
286      !
287      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
288         DO jj = 2, jpjm1   
289            DO ji = fs_2, fs_jpim1   ! vector opt.
290               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
291            END DO
292         END DO
293      END DO
294      !
295      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
296         DO ji = fs_2, fs_jpim1   ! vector opt.
297            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
298            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
299               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
300         END DO
301      END DO
302      DO jk = 2, jpkm1
303         DO jj = 2, jpjm1
304            DO ji = fs_2, fs_jpim1
305               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
306            END DO
307         END DO
308      END DO
309      !
310      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
311         DO ji = fs_2, fs_jpim1   ! vector opt.
312            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
313         END DO
314      END DO
315      DO jk = jpk-2, 1, -1
316         DO jj = 2, jpjm1
317            DO ji = fs_2, fs_jpim1
318               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
319            END DO
320         END DO
321      END DO
322      !
323      !              !==  Vertical diffusion on v  ==!
324      !
325      !                       !* Matrix construction
326      zdt = r2dt * 0.5
327      IF( ln_zad_Aimp ) THEN   !!
328         SELECT CASE( nldf_dyn )
329         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
330            DO jk = 1, jpkm1
331               DO jj = 2, jpjm1 
332                  DO ji = fs_2, fs_jpim1   ! vector opt.
333                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
334                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
335                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
336                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
337                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
338                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
339                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
340                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
341                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
342                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
343                  END DO
344               END DO
345            END DO
346         CASE DEFAULT               ! iso-level lateral mixing
347            DO jk = 1, jpkm1
348               DO jj = 2, jpjm1 
349                  DO ji = fs_2, fs_jpim1   ! vector opt.
350                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
351                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
352                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
353                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
354                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
355                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
356                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
357                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
358                  END DO
359               END DO
360            END DO
361         END SELECT
362         DO jj = 2, jpjm1     !* Surface boundary conditions
363            DO ji = fs_2, fs_jpim1   ! vector opt.
364               zwi(ji,jj,1) = 0._wp
365               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
366               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
367               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) )
368               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
369               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
370            END DO
371         END DO
372      ELSE
373         SELECT CASE( nldf_dyn )
374         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
375            DO jk = 1, jpkm1
376               DO jj = 2, jpjm1   
377                  DO ji = fs_2, fs_jpim1   ! vector opt.
378                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
379                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
380                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
381                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
382                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
383                     zwi(ji,jj,jk) = zzwi
384                     zws(ji,jj,jk) = zzws
385                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
386                  END DO
387               END DO
388            END DO
389         CASE DEFAULT               ! iso-level lateral mixing
390            DO jk = 1, jpkm1
391               DO jj = 2, jpjm1   
392                  DO ji = fs_2, fs_jpim1   ! vector opt.
393                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
394                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
395                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
396                     zwi(ji,jj,jk) = zzwi
397                     zws(ji,jj,jk) = zzws
398                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
399                  END DO
400               END DO
401            END DO
402         END SELECT
403         DO jj = 2, jpjm1        !* Surface boundary conditions
404            DO ji = fs_2, fs_jpim1   ! vector opt.
405               zwi(ji,jj,1) = 0._wp
406               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
407            END DO
408         END DO
409      ENDIF
410      !
411      !              !==  Apply semi-implicit top/bottom friction  ==!
412      !
413      !     Only needed for semi-implicit bottom friction setup. The explicit
414      !     bottom friction has been included in "u(v)a" which act as the R.H.S
415      !     column vector of the tri-diagonal matrix equation
416      !
417      IF( ln_drgimp ) THEN
418         DO jj = 2, jpjm1
419            DO ji = 2, jpim1
420               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
421               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
422               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
423            END DO
424         END DO
425         IF ( ln_isfcav ) THEN
426            DO jj = 2, jpjm1
427               DO ji = 2, jpim1
428                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
429                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
430                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
431               END DO
432            END DO
433         ENDIF
434      ENDIF
435
436      ! Matrix inversion
437      !-----------------------------------------------------------------------
438      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
439      !
440      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
441      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
442      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
443      !        (        ...               )( ...  ) ( ...  )
444      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
445      !
446      !   m is decomposed in the product of an upper and lower triangular matrix
447      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
448      !   The solution (after velocity) is in 2d array va
449      !-----------------------------------------------------------------------
450      !
451      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
452         DO jj = 2, jpjm1   
453            DO ji = fs_2, fs_jpim1   ! vector opt.
454               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
455            END DO
456         END DO
457      END DO
458      !
459      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
460         DO ji = fs_2, fs_jpim1   ! vector opt.         
461            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
462            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
463               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
464         END DO
465      END DO
466      DO jk = 2, jpkm1
467         DO jj = 2, jpjm1
468            DO ji = fs_2, fs_jpim1   ! vector opt.
469               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
470            END DO
471         END DO
472      END DO
473      !
474      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
475         DO ji = fs_2, fs_jpim1   ! vector opt.
476            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
477         END DO
478      END DO
479      DO jk = jpk-2, 1, -1
480         DO jj = 2, jpjm1
481            DO ji = fs_2, fs_jpim1
482               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
483            END DO
484         END DO
485      END DO
486      !
487      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
488         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
489         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
490         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
491         DEALLOCATE( ztrdu, ztrdv ) 
492      ENDIF
493      !                                          ! print mean trends (used for debugging)
494      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
495         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
496         !
497      IF( ln_timing )   CALL timing_stop('dyn_zdf')
498      !
499   END SUBROUTINE dyn_zdf
500
501   !!==============================================================================
502END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.