source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90 @ 9976

Last change on this file since 9976 was 9976, checked in by acc, 3 years ago

Branch: dev_r9956_ENHANCE05_ZAD_AIMP. First set of changes to implement an adaptive-implicit vertical advection option (see ticket #2042). This code compiles and runs but has issues when the new option is activated.

  • Property svn:keywords set to Id
File size: 30.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 licence (./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_e3un, z1_e3vn   !   -      -
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         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==!
163            SELECT CASE( nldf_dyn )
164            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
165               DO jk = 1, jpkm1
166                  DO jj = 2, jpjm1 
167                     DO ji = fs_2, fs_jpim1   ! vector opt.
168                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
169                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm   (ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
170                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
171                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm   (ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
172                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
173                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
174                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
175                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
176                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
177                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
178                     END DO
179                  END DO
180               END DO
181            CASE DEFAULT               ! iso-level lateral mixing
182               DO jk = 1, jpkm1
183                  DO jj = 2, jpjm1 
184                     DO ji = fs_2, fs_jpim1   ! vector opt.
185                        z1_e3un = 1._wp / e3u_n(ji,jj,jk)
186                        zzwi = ( ( avm   (ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
187                           &     / e3uw_a(ji  ,jj,jk  ) ) * z1_e3un * wumask(ji,jj,jk  )
188                        zzws = ( ( avm   (ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
189                           &     / e3uw_a(ji  ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
190                        zWu  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji+1,jj,jk  )   &
191                           &               + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1)   )
192                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
193                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
194                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
195                     END DO
196                  END DO
197               END DO
198            END SELECT
199         ELSE                          !==  Flux form advection  ==!
200            SELECT CASE( nldf_dyn )
201            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
202               DO jk = 1, jpkm1
203                  DO jj = 2, jpjm1 
204                     DO ji = fs_2, fs_jpim1   ! vector opt.
205                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
206                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
207                           &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
208                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
209                           &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
210                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
211                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
212                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
213                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
214                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) )
215                     END DO
216                  END DO
217               END DO
218            CASE DEFAULT               ! iso-level lateral mixing
219               DO jk = 1, jpkm1
220                  DO jj = 2, jpjm1 
221                     DO ji = fs_2, fs_jpim1   ! vector opt.
222                        ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
223                        zzwi = ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
224                        zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
225                        zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
226                        zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
227                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWui, 0._wp ) ) 
228                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWus, 0._wp ) )
229                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) )
230                     END DO
231                  END DO
232               END DO
233            END SELECT
234         ENDIF
235      ELSE
236         SELECT CASE( nldf_dyn )
237         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
238            DO jk = 1, jpkm1
239               DO jj = 2, jpjm1 
240                  DO ji = fs_2, fs_jpim1   ! vector opt.
241                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
242                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
243                        &         / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
244                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
245                        &         / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
246                     zwi(ji,jj,jk) = zzwi
247                     zws(ji,jj,jk) = zzws
248                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
249                  END DO
250               END DO
251            END DO
252         CASE DEFAULT               ! iso-level lateral mixing
253            DO jk = 1, jpkm1
254               DO jj = 2, jpjm1 
255                  DO ji = fs_2, fs_jpim1   ! vector opt.
256                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
257                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_a(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
258                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
259                     zwi(ji,jj,jk) = zzwi
260                     zws(ji,jj,jk) = zzws
261                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
262                  END DO
263               END DO
264            END DO
265         END SELECT
266      ENDIF
267      !
268      DO jj = 2, jpjm1     !* Surface boundary conditions
269         DO ji = fs_2, fs_jpim1   ! vector opt.
270            zwi(ji,jj,1) = 0._wp
271            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
272         END DO
273      END DO
274      !
275      !              !==  Apply semi-implicit bottom friction  ==!
276      !
277      !     Only needed for semi-implicit bottom friction setup. The explicit
278      !     bottom friction has been included in "u(v)a" which act as the R.H.S
279      !     column vector of the tri-diagonal matrix equation
280      !
281      IF ( ln_drgimp ) THEN      ! implicit bottom friction
282         DO jj = 2, jpjm1
283            DO ji = 2, jpim1
284               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
285               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
286               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
287            END DO
288         END DO
289         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
290            DO jj = 2, jpjm1
291               DO ji = 2, jpim1
292                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
293                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
294                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
295                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
296               END DO
297            END DO
298         END IF
299      ENDIF
300      !
301      ! Matrix inversion starting from the first level
302      !-----------------------------------------------------------------------
303      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
304      !
305      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
306      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
307      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
308      !        (        ...               )( ...  ) ( ...  )
309      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
310      !
311      !   m is decomposed in the product of an upper and a lower triangular matrix
312      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
313      !   The solution (the after velocity) is in ua
314      !-----------------------------------------------------------------------
315      !
316      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
317         DO jj = 2, jpjm1   
318            DO ji = fs_2, fs_jpim1   ! vector opt.
319               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
320            END DO
321         END DO
322      END DO
323      !
324      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
325         DO ji = fs_2, fs_jpim1   ! vector opt.
326            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
327            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
328               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
329         END DO
330      END DO
331      DO jk = 2, jpkm1
332         DO jj = 2, jpjm1
333            DO ji = fs_2, fs_jpim1
334               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
335            END DO
336         END DO
337      END DO
338      !
339      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
340         DO ji = fs_2, fs_jpim1   ! vector opt.
341            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
342         END DO
343      END DO
344      DO jk = jpk-2, 1, -1
345         DO jj = 2, jpjm1
346            DO ji = fs_2, fs_jpim1
347               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
348            END DO
349         END DO
350      END DO
351      !
352      !              !==  Vertical diffusion on v  ==!
353      !
354      !                       !* Matrix construction
355      zdt = r2dt * 0.5
356      IF( ln_zad_Aimp ) THEN   !!
357         IF( ln_dynadv_vec ) THEN      !==  Vector invariant advection  ==!
358            SELECT CASE( nldf_dyn )
359            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
360               DO jk = 1, jpkm1
361                  DO jj = 2, jpjm1 
362                     DO ji = fs_2, fs_jpim1   ! vector opt.
363                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
364                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm   (ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
365                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
366                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm   (ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
367                           &     / e3vw_a(ji,jj  ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
368                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
369                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
370                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
371                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
372                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
373                     END DO
374                  END DO
375               END DO
376            CASE DEFAULT               ! iso-level lateral mixing
377               DO jk = 1, jpkm1
378                  DO jj = 2, jpjm1 
379                     DO ji = fs_2, fs_jpim1   ! vector opt.
380                        z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
381                        zzwi = ( ( avm   (ji,jj+1,jk  ) + avm(ji,jj,jk  ) )   &
382                           &     / e3vw_a(ji,jj  ,jk  ) ) * z1_e3vn * wvmask(ji,jj,jk  )
383                        zzws = ( ( avm   (ji,jj+1,jk+1) + avm(ji,jj,jk+1) )   &
384                           &     / e3vw_a(ji  ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
385                        zWv  = 0.25_wp * (   wi(ji,jj,jk  ) + wi(ji,jj+1,jk  )   &
386                           &               + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1)   )
387                        zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
388                        zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
389                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
390                     END DO
391                  END DO
392               END DO
393            END SELECT
394         ELSE                          !==  Flux form advection  ==!
395            SELECT CASE( nldf_dyn )
396            CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
397               DO jk = 1, jpkm1
398                  DO jj = 2, jpjm1 
399                     DO ji = fs_2, fs_jpim1   ! vector opt.
400                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point
401                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
402                           &         / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
403                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
404                           &         / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
405                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
406                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
407                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) ) 
408                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) )
409                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
410                     END DO
411                  END DO
412               END DO
413            CASE DEFAULT               ! iso-level lateral mixing
414               DO jk = 1, jpkm1
415                  DO jj = 2, jpjm1 
416                     DO ji = fs_2, fs_jpim1   ! vector opt.
417                        ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at U-point
418                        zzwi = ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_a(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
419                        zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
420                        zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) )
421                        zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
422                        zwi(ji,jj,jk) = - zdt * ( zzwi +  MIN( zWvi, 0._wp ) ) 
423                        zws(ji,jj,jk) = - zdt * ( zzws -  MAX( zWvs, 0._wp ) )
424                        zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
425                     END DO
426                  END DO
427               END DO
428            END SELECT
429         ENDIF
430      ELSE
431         SELECT CASE( nldf_dyn )
432         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
433            DO jk = 1, jpkm1
434               DO jj = 2, jpjm1   
435                  DO ji = fs_2, fs_jpim1   ! vector opt.
436                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
437                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
438                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
439                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
440                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
441                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
442                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
443                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
444                  END DO
445               END DO
446            END DO
447         CASE DEFAULT               ! iso-level lateral mixing
448            DO jk = 1, jpkm1
449               DO jj = 2, jpjm1   
450                  DO ji = fs_2, fs_jpim1   ! vector opt.
451                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
452                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
453                     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)
454                     zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
455                     zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
456                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
457                  END DO
458               END DO
459            END DO
460         END SELECT
461      ENDIF
462      !
463      DO jj = 2, jpjm1        !* Surface boundary conditions
464         DO ji = fs_2, fs_jpim1   ! vector opt.
465            zwi(ji,jj,1) = 0._wp
466            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
467         END DO
468      END DO
469      !              !==  Apply semi-implicit top/bottom friction  ==!
470      !
471      !     Only needed for semi-implicit bottom friction setup. The explicit
472      !     bottom friction has been included in "u(v)a" which act as the R.H.S
473      !     column vector of the tri-diagonal matrix equation
474      !
475      IF( ln_drgimp ) THEN
476         DO jj = 2, jpjm1
477            DO ji = 2, jpim1
478               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
479               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
480               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
481            END DO
482         END DO
483         IF ( ln_isfcav ) THEN
484            DO jj = 2, jpjm1
485               DO ji = 2, jpim1
486                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
487                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
488                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
489               END DO
490            END DO
491         ENDIF
492      ENDIF
493
494      ! Matrix inversion
495      !-----------------------------------------------------------------------
496      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
497      !
498      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
499      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
500      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
501      !        (        ...               )( ...  ) ( ...  )
502      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
503      !
504      !   m is decomposed in the product of an upper and lower triangular matrix
505      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
506      !   The solution (after velocity) is in 2d array va
507      !-----------------------------------------------------------------------
508      !
509      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
510         DO jj = 2, jpjm1   
511            DO ji = fs_2, fs_jpim1   ! vector opt.
512               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
513            END DO
514         END DO
515      END DO
516      !
517      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
518         DO ji = fs_2, fs_jpim1   ! vector opt.         
519            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
520            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
521               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
522         END DO
523      END DO
524      DO jk = 2, jpkm1
525         DO jj = 2, jpjm1
526            DO ji = fs_2, fs_jpim1   ! vector opt.
527               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
528            END DO
529         END DO
530      END DO
531      !
532      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
533         DO ji = fs_2, fs_jpim1   ! vector opt.
534            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
535         END DO
536      END DO
537      DO jk = jpk-2, 1, -1
538         DO jj = 2, jpjm1
539            DO ji = fs_2, fs_jpim1
540               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
541            END DO
542         END DO
543      END DO
544      !
545      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
546         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
547         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
548         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
549         DEALLOCATE( ztrdu, ztrdv ) 
550      ENDIF
551      !                                          ! print mean trends (used for debugging)
552      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
553         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
554         !
555      IF( ln_timing )   CALL timing_stop('dyn_zdf')
556      !
557   END SUBROUTINE dyn_zdf
558
559   !!==============================================================================
560END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.