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

source: NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90 @ 13470

Last change on this file since 13470 was 13470, checked in by smasson, 4 years ago

r4_trunk: second change of DO loops for routines to be merged, see #2523
  • Property svn:keywords set to Id
File size: 22.6 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) = ( uu(:,:,jk,Nnn) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
113            va(:,:,jk) = ( vv(:,:,jk,Nnn) + 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) * uu(:,:,jk,Nnn)  &
118               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
119            va(:,:,jk) = (         e3v_b(:,:,jk) * vv(:,:,jk,Nnn)  &
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_2D( 0, 0, 0, 0 )
134            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
135            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
136            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
137            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
138            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
139            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
140         END_2D
141         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
142            DO_2D( 0, 0, 0, 0 )
143               iku = miku(ji,jj)         ! top ocean level at u- and v-points
144               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
145               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
146               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
147               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
148               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
149            END_2D
150         END IF
151      ENDIF
152      !
153      !              !==  Vertical diffusion on u  ==!
154      !
155      !                    !* Matrix construction
156      zdt = r2dt * 0.5
157      IF( ln_zad_Aimp ) THEN   !!
158         SELECT CASE( nldf_dyn )
159         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
160            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
161               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
162               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
163                  &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
164               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
165                  &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
166               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
167               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
168               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
169               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
170               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
171            END_3D
172         CASE DEFAULT               ! iso-level lateral mixing
173            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
174               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
175               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
176               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)
177               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
178               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
179               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
180               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
181               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
182            END_3D
183         END SELECT
184         DO_2D( 0, 0, 0, 0 )
185            zwi(ji,jj,1) = 0._wp
186            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
187            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
188            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
189            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
190            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
191         END_2D
192      ELSE
193         SELECT CASE( nldf_dyn )
194         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
195            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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  ) + akzu(ji,jj,jk  ) )   &
198                  &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
199               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
200                  &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
201               zwi(ji,jj,jk) = zzwi
202               zws(ji,jj,jk) = zzws
203               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
204            END_3D
205         CASE DEFAULT               ! iso-level lateral mixing
206            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
207               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at U-point
208               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
209               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)
210               zwi(ji,jj,jk) = zzwi
211               zws(ji,jj,jk) = zzws
212               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
213            END_3D
214         END SELECT
215         DO_2D( 0, 0, 0, 0 )
216            zwi(ji,jj,1) = 0._wp
217            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
218         END_2D
219      ENDIF
220      !
221      !
222      !              !==  Apply semi-implicit bottom friction  ==!
223      !
224      !     Only needed for semi-implicit bottom friction setup. The explicit
225      !     bottom friction has been included in "u(v)a" which act as the R.H.S
226      !     column vector of the tri-diagonal matrix equation
227      !
228      IF ( ln_drgimp ) THEN      ! implicit bottom friction
229         DO_2D( 0, 0, 0, 0 )
230            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
231            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
232            zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
233         END_2D
234         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
235            DO_2D( 0, 0, 0, 0 )
236               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
237               iku = miku(ji,jj)       ! ocean top level at u- and v-points
238               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
239               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
240            END_2D
241         END IF
242      ENDIF
243      !
244      ! Matrix inversion starting from the first level
245      !-----------------------------------------------------------------------
246      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
247      !
248      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
249      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
250      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
251      !        (        ...               )( ...  ) ( ...  )
252      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
253      !
254      !   m is decomposed in the product of an upper and a lower triangular matrix
255      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
256      !   The solution (the after velocity) is in ua
257      !-----------------------------------------------------------------------
258      !
259      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
260         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
261      END_3D
262      !
263      DO_2D( 0, 0, 0, 0 )
264         ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
265         ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
266            &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
267      END_2D
268      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
269         ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
270      END_3D
271      !
272      DO_2D( 0, 0, 0, 0 )
273         ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
274      END_2D
275      DO_3D( 0, 0, 0, 0, jpk-2, 1, -1 )
276         ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
277      END_3D
278      !
279      !              !==  Vertical diffusion on v  ==!
280      !
281      !                       !* Matrix construction
282      zdt = r2dt * 0.5
283      IF( ln_zad_Aimp ) THEN   !!
284         SELECT CASE( nldf_dyn )
285         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
286            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
287               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
288               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
289                  &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
290               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
291                  &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
292               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
293               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
294               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
295               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
296               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
297            END_3D
298         CASE DEFAULT               ! iso-level lateral mixing
299            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
300               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
301               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
302               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)
303               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
304               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
305               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
306               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
307               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
308            END_3D
309         END SELECT
310         DO_2D( 0, 0, 0, 0 )
311            zwi(ji,jj,1) = 0._wp
312            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
313            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
314            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
315            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
316            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
317         END_2D
318      ELSE
319         SELECT CASE( nldf_dyn )
320         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
321            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
322               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at V-point
323               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
324                  &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
325               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
326                  &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
327               zwi(ji,jj,jk) = zzwi
328               zws(ji,jj,jk) = zzws
329               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
330            END_3D
331         CASE DEFAULT               ! iso-level lateral mixing
332            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
335               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)
336               zwi(ji,jj,jk) = zzwi
337               zws(ji,jj,jk) = zzws
338               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
339            END_3D
340         END SELECT
341         DO_2D( 0, 0, 0, 0 )
342            zwi(ji,jj,1) = 0._wp
343            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
344         END_2D
345      ENDIF
346      !
347      !              !==  Apply semi-implicit top/bottom friction  ==!
348      !
349      !     Only needed for semi-implicit bottom friction setup. The explicit
350      !     bottom friction has been included in "u(v)a" which act as the R.H.S
351      !     column vector of the tri-diagonal matrix equation
352      !
353      IF( ln_drgimp ) THEN
354         DO_2D( 0, 0, 0, 0 )
355            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
356            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
357            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
358         END_2D
359         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
360            DO_2D( 0, 0, 0, 0 )
361               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
362               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
363               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
364            END_2D
365         ENDIF
366      ENDIF
367
368      ! Matrix inversion
369      !-----------------------------------------------------------------------
370      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
371      !
372      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
373      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
374      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
375      !        (        ...               )( ...  ) ( ...  )
376      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
377      !
378      !   m is decomposed in the product of an upper and lower triangular matrix
379      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
380      !   The solution (after velocity) is in 2d array va
381      !-----------------------------------------------------------------------
382      !
383      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
384         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
385      END_3D
386      !
387      DO_2D( 0, 0, 0, 0 )
388         ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
389         va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
390            &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
391      END_2D
392      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
393         va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
394      END_3D
395      !
396      DO_2D( 0, 0, 0, 0 )
397         va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
398      END_2D
399      DO_3D( 0, 0, 0, 0, jpk-2, 1, -1 )
400         va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
401      END_3D
402      !
403      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
404         ztrdu(:,:,:) = ( ua(:,:,:) - uu(:,:,:,Nnn) ) / r2dt - ztrdu(:,:,:)
405         ztrdv(:,:,:) = ( va(:,:,:) - vv(:,:,:,Nnn) ) / r2dt - ztrdv(:,:,:)
406         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
407         DEALLOCATE( ztrdu, ztrdv ) 
408      ENDIF
409      !                                          ! print mean trends (used for debugging)
410      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
411         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
412         !
413      IF( ln_timing )   CALL timing_stop('dyn_zdf')
414      !
415   END SUBROUTINE dyn_zdf
416
417   !!==============================================================================
418END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.