source: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynzdf.F90 @ 11613

Last change on this file since 11613 was 11613, checked in by davestorkey, 14 months ago

UKMO/NEMO_4.0_momentum_trends branch : first set of code changes.

File size: 27.5 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 iom             ! IOM library
28   USE lib_mpp        ! MPP library
29   USE prtctl         ! Print control
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_zdf   !  routine called by step.F90
36
37   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
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 license (see ./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      INTEGER  ::   ikbu, ikbv         ! local integers
73      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
74      REAL(wp) ::   zzws, ze3va        !   -      -
75      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
76      REAL(wp) ::   zWu , zWv          !   -      -
77      REAL(wp) ::   zWui, zWvi         !   -      -
78      REAL(wp) ::   zWus, zWvs         !   -      -
79      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace
80      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
81      !!---------------------------------------------------------------------
82      !
83      IF( ln_timing )   CALL timing_start('dyn_zdf')
84      !
85      IF( kt == nit000 ) THEN       !* initialization
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
89         !
90         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
91         ELSE                   ;    r_vvl = 1._wp
92         ENDIF
93      ENDIF
94      !                             !* set time step
95      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
96      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
97      ENDIF
98      !
99      !                             !* explicit top/bottom drag case
100      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, ub, vb, ua, va )  ! add top/bottom friction trend to (ua,va)
101      !
102      !
103      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
104      !
105      !                    ! time stepping except vertical diffusion
106      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
107         DO jk = 1, jpkm1
108            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
109            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
110         END DO
111      ELSE                                      ! applied on thickness weighted velocity
112         DO jk = 1, jpkm1
113            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
114               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
115            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
116               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
117         END DO
118      ENDIF
119      !                    ! add top/bottom friction
120      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
121      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
122      !                not lead to the effective stress seen over the whole barotropic loop.
123      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
124      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
125         DO jk = 1, jpkm1        ! remove barotropic velocities
126            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
127            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
128         END DO
129         IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
130            ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
131            ztrdu(:,:,:) = ua(:,:,:)
132            ztrdv(:,:,:) = va(:,:,:)
133         ENDIF
134         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
137               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
138               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
139               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
140               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
141               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
142            END DO
143         END DO
144         IF( l_trddyn )   THEN                      ! save explicit part of bottom friction trends
145            ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt 
146            ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt 
147            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfr, kt )
148         ENDIF
149         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
150            DO jj = 2, jpjm1       
151               DO ji = fs_2, fs_jpim1   ! vector opt.
152                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
153                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
154                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
155                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
156                  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
157                  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
158               END DO
159            END DO
160         END IF
161      ENDIF
162      !
163      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
164         IF( .NOT. ALLOCATED(ztrdu) ) ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
165         ztrdu(:,:,:) = ua(:,:,:)
166         ztrdv(:,:,:) = va(:,:,:)
167      ENDIF
168      !
169      !              !==  Vertical diffusion on u  ==!
170      !
171      !                    !* Matrix construction
172      zdt = r2dt * 0.5
173      IF( ln_zad_Aimp ) THEN   !!
174         SELECT CASE( nldf_dyn )
175         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
176            DO jk = 1, jpkm1
177               DO jj = 2, jpjm1 
178                  DO ji = fs_2, fs_jpim1   ! vector opt.
179                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
180                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
181                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
182                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
183                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
184                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
185                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
186                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
187                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
188                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
189                  END DO
190               END DO
191            END DO
192         CASE DEFAULT               ! iso-level lateral mixing
193            DO jk = 1, jpkm1
194               DO jj = 2, jpjm1 
195                  DO ji = fs_2, fs_jpim1   ! vector opt.
196                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
197                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
198                     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)
199                     zWui = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) )
200                     zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
201                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
202                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
203                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
204                  END DO
205               END DO
206            END DO
207         END SELECT
208         DO jj = 2, jpjm1     !* Surface boundary conditions
209            DO ji = fs_2, fs_jpim1   ! vector opt.
210               zwi(ji,jj,1) = 0._wp
211               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
212               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
213               zWus = 0.5_wp * ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) )
214               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
215               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
216            END DO
217         END DO
218      ELSE
219         SELECT CASE( nldf_dyn )
220         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
221            DO jk = 1, jpkm1
222               DO jj = 2, jpjm1 
223                  DO ji = fs_2, fs_jpim1   ! vector opt.
224                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
225                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
226                        &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
227                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
228                        &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
229                     zwi(ji,jj,jk) = zzwi
230                     zws(ji,jj,jk) = zzws
231                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
232                  END DO
233               END DO
234            END DO
235         CASE DEFAULT               ! iso-level lateral mixing
236            DO jk = 1, jpkm1
237               DO jj = 2, jpjm1 
238                  DO ji = fs_2, fs_jpim1   ! vector opt.
239                     ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
240                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
241                     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)
242                     zwi(ji,jj,jk) = zzwi
243                     zws(ji,jj,jk) = zzws
244                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
245                  END DO
246               END DO
247            END DO
248         END SELECT
249         DO jj = 2, jpjm1     !* Surface boundary conditions
250            DO ji = fs_2, fs_jpim1   ! vector opt.
251               zwi(ji,jj,1) = 0._wp
252               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
253            END DO
254         END DO
255      ENDIF
256      !
257      !
258      !              !==  Apply semi-implicit bottom friction  ==!
259      !
260      !     Only needed for semi-implicit bottom friction setup. The explicit
261      !     bottom friction has been included in "u(v)a" which act as the R.H.S
262      !     column vector of the tri-diagonal matrix equation
263      !
264      IF ( ln_drgimp ) THEN      ! implicit bottom friction
265         DO jj = 2, jpjm1
266            DO ji = 2, jpim1
267               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
268               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
269               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
270            END DO
271         END DO
272         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
273            DO jj = 2, jpjm1
274               DO ji = 2, jpim1
275                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
276                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
277                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
278                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
279               END DO
280            END DO
281         END IF
282      ENDIF
283      !
284      ! Matrix inversion starting from the first level
285      !-----------------------------------------------------------------------
286      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
287      !
288      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
289      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
290      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
291      !        (        ...               )( ...  ) ( ...  )
292      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
293      !
294      !   m is decomposed in the product of an upper and a lower triangular matrix
295      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
296      !   The solution (the after velocity) is in ua
297      !-----------------------------------------------------------------------
298      !
299      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
300         DO jj = 2, jpjm1   
301            DO ji = fs_2, fs_jpim1   ! vector opt.
302               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
303            END DO
304         END DO
305      END DO
306      !
307      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
308         DO ji = fs_2, fs_jpim1   ! vector opt.
309            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
310            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
311               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
312         END DO
313      END DO
314      DO jk = 2, jpkm1
315         DO jj = 2, jpjm1
316            DO ji = fs_2, fs_jpim1
317               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
318            END DO
319         END DO
320      END DO
321      !
322      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
323         DO ji = fs_2, fs_jpim1   ! vector opt.
324            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
325         END DO
326      END DO
327      DO jk = jpk-2, 1, -1
328         DO jj = 2, jpjm1
329            DO ji = fs_2, fs_jpim1
330               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
331            END DO
332         END DO
333      END DO
334      !
335      !              !==  Vertical diffusion on v  ==!
336      !
337      !                       !* Matrix construction
338      zdt = r2dt * 0.5
339      IF( ln_zad_Aimp ) THEN   !!
340         SELECT CASE( nldf_dyn )
341         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
342            DO jk = 1, jpkm1
343               DO jj = 2, jpjm1 
344                  DO ji = fs_2, fs_jpim1   ! vector opt.
345                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
346                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
347                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
348                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
349                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
350                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
351                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
352                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
353                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
354                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
355                  END DO
356               END DO
357            END DO
358         CASE DEFAULT               ! iso-level lateral mixing
359            DO jk = 1, jpkm1
360               DO jj = 2, jpjm1 
361                  DO ji = fs_2, fs_jpim1   ! vector opt.
362                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
363                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
364                     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)
365                     zWvi = 0.5_wp * ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) * wvmask(ji,jj,jk  )
366                     zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
367                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
368                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
369                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
370                  END DO
371               END DO
372            END DO
373         END SELECT
374         DO jj = 2, jpjm1     !* Surface boundary conditions
375            DO ji = fs_2, fs_jpim1   ! vector opt.
376               zwi(ji,jj,1) = 0._wp
377               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
378               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
379               zWvs = 0.5_wp * ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) )
380               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
381               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
382            END DO
383         END DO
384      ELSE
385         SELECT CASE( nldf_dyn )
386         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
387            DO jk = 1, jpkm1
388               DO jj = 2, jpjm1   
389                  DO ji = fs_2, fs_jpim1   ! vector opt.
390                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
391                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
392                        &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
393                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
394                        &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
395                     zwi(ji,jj,jk) = zzwi
396                     zws(ji,jj,jk) = zzws
397                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
398                  END DO
399               END DO
400            END DO
401         CASE DEFAULT               ! iso-level lateral mixing
402            DO jk = 1, jpkm1
403               DO jj = 2, jpjm1   
404                  DO ji = fs_2, fs_jpim1   ! vector opt.
405                     ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
406                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
407                     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)
408                     zwi(ji,jj,jk) = zzwi
409                     zws(ji,jj,jk) = zzws
410                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
411                  END DO
412               END DO
413            END DO
414         END SELECT
415         DO jj = 2, jpjm1        !* Surface boundary conditions
416            DO ji = fs_2, fs_jpim1   ! vector opt.
417               zwi(ji,jj,1) = 0._wp
418               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
419            END DO
420         END DO
421      ENDIF
422      !
423      !              !==  Apply semi-implicit top/bottom friction  ==!
424      !
425      !     Only needed for semi-implicit bottom friction setup. The explicit
426      !     bottom friction has been included in "u(v)a" which act as the R.H.S
427      !     column vector of the tri-diagonal matrix equation
428      !
429      IF( ln_drgimp ) THEN
430         DO jj = 2, jpjm1
431            DO ji = 2, jpim1
432               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
433               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
434               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
435            END DO
436         END DO
437         IF ( ln_isfcav ) THEN
438            DO jj = 2, jpjm1
439               DO ji = 2, jpim1
440                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
441                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
442                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
443               END DO
444            END DO
445         ENDIF
446      ENDIF
447
448      ! Matrix inversion
449      !-----------------------------------------------------------------------
450      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
451      !
452      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
453      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
454      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
455      !        (        ...               )( ...  ) ( ...  )
456      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
457      !
458      !   m is decomposed in the product of an upper and lower triangular matrix
459      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
460      !   The solution (after velocity) is in 2d array va
461      !-----------------------------------------------------------------------
462      !
463      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
464         DO jj = 2, jpjm1   
465            DO ji = fs_2, fs_jpim1   ! vector opt.
466               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
467            END DO
468         END DO
469      END DO
470      !
471      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
472         DO ji = fs_2, fs_jpim1   ! vector opt.         
473            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
474            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
475               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
476         END DO
477      END DO
478      DO jk = 2, jpkm1
479         DO jj = 2, jpjm1
480            DO ji = fs_2, fs_jpim1   ! vector opt.
481               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
482            END DO
483         END DO
484      END DO
485      !
486      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
487         DO ji = fs_2, fs_jpim1   ! vector opt.
488            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
489         END DO
490      END DO
491      DO jk = jpk-2, 1, -1
492         DO jj = 2, jpjm1
493            DO ji = fs_2, fs_jpim1
494               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
495            END DO
496         END DO
497      END DO
498      !
499      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
500
501         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt 
502         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt 
503         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
504         !
505         IF( ln_drgimp ) THEN
506            ztrdu(:,:,:) = 0._wp    ;   ztrdv(:,:,:) = 0._wp
507            DO jk = 1, jpkm1
508               DO jj = 2, jpjm1
509                  DO ji = 2, jpim1
510                     ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
511                     ikbv = mbkv(ji,jj)
512                     ztrdu(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) )  & 
513     &                               * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
514                     ztrdv(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) )  &
515     &                               * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
516                  END DO
517               END DO
518            END DO
519            CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfri, kt )
520         ENDIF
521         !
522         DEALLOCATE( ztrdu, ztrdv ) 
523      ENDIF
524      !                                          ! print mean trends (used for debugging)
525      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
526         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
527         !
528      IF( ln_timing )   CALL timing_stop('dyn_zdf')
529      !
530   END SUBROUTINE dyn_zdf
531
532   !!==============================================================================
533END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.