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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynzdf.F90 @ 12293

Last change on this file since 12293 was 12293, checked in by mathiot, 4 years ago

part of ticket #2358 also affect dev_r11943_MERGE_2019 (option2 seems OK)
  • Property svn:keywords set to Id
File size: 27.0 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, Kbb, Kmm, Krhs, puu, pvv, Kaa )
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      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
57      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after)   otherwise
58      !!               - update the after velocity with the implicit vertical mixing.
59      !!      This requires to solver the following system:
60      !!         u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) 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 :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity
66      !!---------------------------------------------------------------------
67      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
68      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
69      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
70      !
71      INTEGER  ::   ji, jj, jk         ! dummy loop indices
72      INTEGER  ::   iku, ikv           ! 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, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
101      !
102      !
103      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
104         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
105         ztrdu(:,:,:) = puu(:,:,:,Krhs)
106         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
107      ENDIF
108      !
109      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
110      !
111      !                    ! time stepping except vertical diffusion
112      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
113         DO jk = 1, jpkm1
114            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)
115            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)
116         END DO
117      ELSE                                      ! applied on thickness weighted velocity
118         DO jk = 1, jpkm1
119            puu(:,:,jk,Kaa) = (         e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb)  &
120               &          + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs)  ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)
121            pvv(:,:,jk,Kaa) = (         e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb)  &
122               &          + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs)  ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)
123         END DO
124      ENDIF
125      !                    ! add top/bottom friction
126      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
127      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
128      !                not lead to the effective stress seen over the whole barotropic loop.
129      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
130      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
131         DO jk = 1, jpkm1        ! remove barotropic velocities
132            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)
133            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)
134         END DO
135         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
136            DO ji = fs_2, fs_jpim1   ! vector opt.
137               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
138               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
139               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
140               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
141               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
142               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
143            END DO
144         END DO
145         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
146            DO jj = 2, jpjm1       
147               DO ji = fs_2, fs_jpim1   ! vector opt.
148                  iku = miku(ji,jj)         ! top ocean level at u- and v-points
149                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
150                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
151                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
152                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
153                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
154               END DO
155            END DO
156         END IF
157      ENDIF
158      !
159      !              !==  Vertical diffusion on u  ==!
160      !
161      !                    !* Matrix construction
162      zdt = r2dt * 0.5
163      IF( ln_zad_Aimp ) THEN   !!
164         SELECT CASE( nldf_dyn )
165         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
166            DO jk = 1, jpkm1
167               DO jj = 2, jpjm1 
168                  DO ji = fs_2, fs_jpim1   ! vector opt.
169                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
170                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
171                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
172                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
173                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
174                     zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
175                     zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
176                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
177                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
178                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
179                  END DO
180               END DO
181            END DO
182         CASE DEFAULT               ! iso-level lateral mixing
183            DO jk = 1, jpkm1
184               DO jj = 2, jpjm1 
185                  DO ji = fs_2, fs_jpim1   ! vector opt.
186                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
187                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
188                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
189                     zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
190                     zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
191                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
192                     zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
193                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
194                  END DO
195               END DO
196            END DO
197         END SELECT
198         DO jj = 2, jpjm1     !* Surface boundary conditions
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               zwi(ji,jj,1) = 0._wp
201               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa)
202               zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
203               zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
204               zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
205               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
206            END DO
207         END DO
208      ELSE
209         SELECT CASE( nldf_dyn )
210         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
211            DO jk = 1, jpkm1
212               DO jj = 2, jpjm1 
213                  DO ji = fs_2, fs_jpim1   ! vector opt.
214                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
215                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
216                        &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
217                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
218                        &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
219                     zwi(ji,jj,jk) = zzwi
220                     zws(ji,jj,jk) = zzws
221                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
222                  END DO
223               END DO
224            END DO
225         CASE DEFAULT               ! iso-level lateral mixing
226            DO jk = 1, jpkm1
227               DO jj = 2, jpjm1 
228                  DO ji = fs_2, fs_jpim1   ! vector opt.
229                     ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
230                     zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
231                     zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
232                     zwi(ji,jj,jk) = zzwi
233                     zws(ji,jj,jk) = zzws
234                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
235                  END DO
236               END DO
237            END DO
238         END SELECT
239         DO jj = 2, jpjm1     !* Surface boundary conditions
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zwi(ji,jj,1) = 0._wp
242               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
243            END DO
244         END DO
245      ENDIF
246      !
247      !
248      !              !==  Apply semi-implicit bottom friction  ==!
249      !
250      !     Only needed for semi-implicit bottom friction setup. The explicit
251      !     bottom friction has been included in "u(v)a" which act as the R.H.S
252      !     column vector of the tri-diagonal matrix equation
253      !
254      IF ( ln_drgimp ) THEN      ! implicit bottom friction
255         DO jj = 2, jpjm1
256            DO ji = 2, jpim1
257               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
258               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
259               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
260            END DO
261         END DO
262         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
263            DO jj = 2, jpjm1
264               DO ji = 2, jpim1
265                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
266                  iku = miku(ji,jj)       ! ocean top level at u- and v-points
267                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
268                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
269               END DO
270            END DO
271         END IF
272      ENDIF
273      !
274      ! Matrix inversion starting from the first level
275      !-----------------------------------------------------------------------
276      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
277      !
278      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
279      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
280      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
281      !        (        ...               )( ...  ) ( ...  )
282      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
283      !
284      !   m is decomposed in the product of an upper and a lower triangular matrix
285      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
286      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
287      !-----------------------------------------------------------------------
288      !
289      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
290         DO jj = 2, jpjm1   
291            DO ji = fs_2, fs_jpim1   ! vector opt.
292               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
293            END DO
294         END DO
295      END DO
296      !
297      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
298         DO ji = fs_2, fs_jpim1   ! vector opt.
299            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 
300            puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
301               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
302         END DO
303      END DO
304      DO jk = 2, jpkm1
305         DO jj = 2, jpjm1
306            DO ji = fs_2, fs_jpim1
307               puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
308            END DO
309         END DO
310      END DO
311      !
312      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
313         DO ji = fs_2, fs_jpim1   ! vector opt.
314            puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
315         END DO
316      END DO
317      DO jk = jpk-2, 1, -1
318         DO jj = 2, jpjm1
319            DO ji = fs_2, fs_jpim1
320               puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
321            END DO
322         END DO
323      END DO
324      !
325      !              !==  Vertical diffusion on v  ==!
326      !
327      !                       !* Matrix construction
328      zdt = r2dt * 0.5
329      IF( ln_zad_Aimp ) THEN   !!
330         SELECT CASE( nldf_dyn )
331         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
332            DO jk = 1, jpkm1
333               DO jj = 2, jpjm1 
334                  DO ji = fs_2, fs_jpim1   ! vector opt.
335                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
336                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
337                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
338                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
339                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
340                     zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
341                     zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
342                     zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
343                     zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
344                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
345                  END DO
346               END DO
347            END DO
348         CASE DEFAULT               ! iso-level lateral mixing
349            DO jk = 1, jpkm1
350               DO jj = 2, jpjm1 
351                  DO ji = fs_2, fs_jpim1   ! vector opt.
352                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
353                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
354                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
355                     zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
356                     zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
357                     zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
358                     zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
359                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
360                  END DO
361               END DO
362            END DO
363         END SELECT
364         DO jj = 2, jpjm1     !* Surface boundary conditions
365            DO ji = fs_2, fs_jpim1   ! vector opt.
366               zwi(ji,jj,1) = 0._wp
367               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa)
368               zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
369               zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
370               zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
371               zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
372            END DO
373         END DO
374      ELSE
375         SELECT CASE( nldf_dyn )
376         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
377            DO jk = 1, jpkm1
378               DO jj = 2, jpjm1   
379                  DO ji = fs_2, fs_jpim1   ! vector opt.
380                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
381                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
382                        &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
383                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
384                        &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
385                     zwi(ji,jj,jk) = zzwi
386                     zws(ji,jj,jk) = zzws
387                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
388                  END DO
389               END DO
390            END DO
391         CASE DEFAULT               ! iso-level lateral mixing
392            DO jk = 1, jpkm1
393               DO jj = 2, jpjm1   
394                  DO ji = fs_2, fs_jpim1   ! vector opt.
395                     ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
396                     zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
397                     zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
398                     zwi(ji,jj,jk) = zzwi
399                     zws(ji,jj,jk) = zzws
400                     zwd(ji,jj,jk) = 1._wp - zzwi - zzws
401                  END DO
402               END DO
403            END DO
404         END SELECT
405         DO jj = 2, jpjm1        !* Surface boundary conditions
406            DO ji = fs_2, fs_jpim1   ! vector opt.
407               zwi(ji,jj,1) = 0._wp
408               zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
409            END DO
410         END DO
411      ENDIF
412      !
413      !              !==  Apply semi-implicit top/bottom friction  ==!
414      !
415      !     Only needed for semi-implicit bottom friction setup. The explicit
416      !     bottom friction has been included in "u(v)a" which act as the R.H.S
417      !     column vector of the tri-diagonal matrix equation
418      !
419      IF( ln_drgimp ) THEN
420         DO jj = 2, jpjm1
421            DO ji = 2, jpim1
422               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
423               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
424               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
425            END DO
426         END DO
427         IF ( ln_isfcav ) THEN
428            DO jj = 2, jpjm1
429               DO ji = 2, jpim1
430                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
431                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
432                  zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
433               END DO
434            END DO
435         ENDIF
436      ENDIF
437
438      ! Matrix inversion
439      !-----------------------------------------------------------------------
440      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
441      !
442      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
443      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
444      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
445      !        (        ...               )( ...  ) ( ...  )
446      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
447      !
448      !   m is decomposed in the product of an upper and lower triangular matrix
449      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
450      !   The solution (after velocity) is in 2d array va
451      !-----------------------------------------------------------------------
452      !
453      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
454         DO jj = 2, jpjm1   
455            DO ji = fs_2, fs_jpim1   ! vector opt.
456               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
457            END DO
458         END DO
459      END DO
460      !
461      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
462         DO ji = fs_2, fs_jpim1   ! vector opt.         
463            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 
464            pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
465               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
466         END DO
467      END DO
468      DO jk = 2, jpkm1
469         DO jj = 2, jpjm1
470            DO ji = fs_2, fs_jpim1   ! vector opt.
471               pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
472            END DO
473         END DO
474      END DO
475      !
476      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
477         DO ji = fs_2, fs_jpim1   ! vector opt.
478            pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
479         END DO
480      END DO
481      DO jk = jpk-2, 1, -1
482         DO jj = 2, jpjm1
483            DO ji = fs_2, fs_jpim1
484               pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
485            END DO
486         END DO
487      END DO
488      !
489      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
490         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:)
491         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:)
492         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
493         DEALLOCATE( ztrdu, ztrdv ) 
494      ENDIF
495      !                                          ! print mean trends (used for debugging)
496      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
497         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
498         !
499      IF( ln_timing )   CALL timing_stop('dyn_zdf')
500      !
501   END SUBROUTINE dyn_zdf
502
503   !!==============================================================================
504END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.